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ABSTRACT 


One of the major concerns in the design of an active control system is obtaining 
the information needed for effective feedback. This involves the combination of 
sensing and estimation. A sensor location index is defined as the weighted sum of 
the mean square estimation errors in which the sensor locations can be regarded 
as estimator design parameters. The design goal is to choose these locations to 
minimize the sensor location index. The gradient of the sensor location index with 
respect to each individual sensor location is formulated and a program using this 
gradient for systematic optimal sensor location search is developed. The choice 
of the number of sensors is a tradeoff between the estimation quality based upon 
the same performance index and the total costs of installing and maintaining extra 
sensors. 

An experimental study for choosing the sensor location is conducted on an 
aeroelastic system. It is the physical realization of a two degree of freedom typical 
section wing. It consists of a NACA 0015 typical section wing with six accelerome- 
ters installed inside along the wing chord as the estimator measuring instruments, an 
existing wind tunnel section, and some other accompanying experimental devices. 
The system modeling which includes the unsteady aerodynamics model developed 
by Stephen Rock has been improved. The center of percussion of the rigid two de- 
gree of freedom typical section wing has been verified as a sensor location for which 
the system is unobservable. Experimental results verify the trend of the theoretical 
predictions of the sensor location index for different sensor locations at various wind 
speeds. 
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Chapter I 


INTRODUCTION 


A. BACKGROUND 

Large structures have been considered for future space missions, such as 
communications, collection of solar energies, etc. Because of the inherent size 
and the necessarily low weight to area ratio, the influence of structural flexibility 
is becoming more significant. This is also true for the future energy efficient air 
transports which in addition need increased performance, comfort, and service 
time. To meet the requirement of larger configurations and more ambitious 
specifications on control system performance, e.g. controlling both the geometric 
shape and the orientation of the configuration (attitude control) [Refs. 1-2], as 
well as vibration suppression [Refs. 3-6], one is led to the use of modern mul- 
tivariable control theory and the associated concepts of controllability and obser- 
vability. However, the concepts of controllability and observability in modern 
control theory can only provide binary information, i.e. either a system is con- 
trollable (observable) or uncontrollable (unobservable) [7]. The more useful 
design questions which are naturally asked still remain unsolved: 

a) What types of actuators and sensors should be used? 

b) How to identify the minimum required number of them, and how does 
the increased number affect the overall performance? 

c) Where should they be placed? 

Among all, placement exhibits a fundamental and largest degree of freedom 
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available to the designer even though it is usually not a very straightforward 
question to answer. 

There are already several papers discussing these problems and some guide- 
lines have been suggested [Refs. 8-14}. 

Juang and Rodriguez [8] start from the general linear quadratic (LQ) for- 
mulations and establish two criteria relating the actuator and sensor placement 
to the minimums of the optimal cost function and state estimation error. This is 
a simple approach but it doesn’t take into account the time constraint of the con- 
trol or estimation. 

Likins, et al. [9] define a measure of controllability based on the minimum 
initial distance from the origin in a state space from which the system can be 
brought to the origin in time T e , what they call “recovery region”. This is a 
difficult problem itself which can only be solved approximately. 

Hughes and Skelton [10] develop the modal controllability (observability) 
by measuring the controllability (observability) norm based on a linear matrix- 
second-order system formulation. 

Kammer and Sesak [11] interpret the increasing insensitivity to parameter 
variations (robustness) versus actuator number analytically, it is analogous to 
that caused by the o-shifted performance index. This result adds to the 
designer’s ability one more degree of freedom, especially in dealing with distri- 
buted systems. It also confirms that the performance tends to increase with 
increasing system complexity, i.e. with more actuators and sensors. 

After reviewing Juang’s and Likins’ works and finding objections, Vander 
Velde and Carignan define the degree of controllability (observability) as a 
linear measure of the weighted “volume” of the hyperellipsoid in the 
transformed equicontrol (equimeasure) space [Refs. 12-14]. This definition will 
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be reviewed in Chapter II. 

It is evident that they all try to find a standardized criterion that measures 
directly the degree of controllability (observability) for different actuator (sen- 
sor) locations. Most important, it must be a quantitative measure that can be 
easily computed and have a direct physical interpretation so that the control sys- 
tem designer can rank the desirability of various candidate actuator (sensor) 
distributions in a meaningful way. 

In this research, our major concern is to find a guideline in this respect for 
sensing and state estimation. So we will hereafter restrict ourselves to the sensor 
part in the above questions. While the duality between estimator and controller 
design is well known and established [Refs. 15-16], this research goal can thus 
easily be extended to the controller part. 

First, a so called sensor location index (J s jr) will be proposed. It is defined 
as a weighted sum of the mean square estimation errors, i.e. J SL = tr ( WP), 
where W is a symmetric positive-definite weighting matrix. It happens to be 
similar to Juang’s idea, except for the weighted sum formulation. The choice of 
W is similar to Vander Velde’s choice of the transformation to an equimeasure 
space, specifically, it is the inverse square of his transformation matrix. 

Using existing software support, it is easy to compute the J SL in a steady 
state manner. Since our goal is to look for the optimal sensor location directly, 
an optimization technique using the gradient search scheme is applied for this 
purpose. The program is designed to report the figure of merit for different sen- 
sor locations and increasing sensor number, thus offers the designer a handy refer- 
ence to make his judgements. 

An experimental study to evaluate the sensor number and their locations is 
conducted on an aeroelastic system. It is the physical realization of a two degree 
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of freedom (DOF) typical section wing [Refs. 17-18]. It consists of a NACA 
0015 typical section wing with six accelerometers installed inside along the wing 
chord as the estimator measuring instruments. The evaluation process is done 
through the Kalman filtering technique for state estimation and sensor location 
index computation. However, effective implementations of Kalman filtering 
depend heavily upon the accuracy of the unsteady aerodynamics theory in 
predicting the loads associated with the general motions of the airlifting surfaces. 
Fortunately, this unsteady aerodynamics theory has already been investigated 
thoroughly. 

A complete solution for the unsteady aerodynamic loads of the typical sec- 
tion wing, including a trailing-edge flap, undergoing a simple harmonic motion 
was first presented by Theodorsen [19] in 1935 for the conditions of incompres- 
sible, invicid, two-dimensional airflows. 

In 1951, Timman [20] solved the same problem except under the con- 
straint of a two-dimensional wind tunnel. 

In 1977, Edwards [Refs. 21-22] extended Theodorsen’s results to include 
arbitrary airfoil motions by recognizing a derivation of the generalized Theodor- 
sen function in the work of Von Karman and Sears [Refs. 23-24]. 

In 1978, Rock [25] combined Timman ’s results with the work of Edwards’ 
to develop an aeroelastic model for a 2 DOF typical section wing undergoing 
arbitrary motions in a small subsonic wind tunnel. He also verified this analytic 
model by experimental investigations. 

Afterwards, in 1981, Stoltz [26] extended Rock’s work by adding a 
trailing-edge flap to the typical section wing. Using this as a control surface, he 
was able to design an active aeroelastic control law capable of stabilizing the 
experimental wing-flap system up to a dynamic pressure equal to twice the open 
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loop flutter value. 

In our experimental studies, Rock’s unsteady aerodynamics modeling is used 
for our new airfoil model and it is presented in Appendix A. 

B. THESIS OUTLINE 

In Chapter II, a continuous free-free beam example is first given to explain 
the physical meaning of system observability. Then the optimal estimator design 
problem is reformulated for correlated process disturbance and measurement 
noise. The derivation is based on minimizing a weighted sum of the mean square 
estimation errors as the performance index [15]. It turns out that it is not trivial 
especially for the discrete case. The OPTSYS program [27] is modified to the 
LOPTSYS to accept a more general input form including a rate measurement. 
It solves for the continuous steady state estimator gain with correlated process 
disturbance and measurement noise. A similar work is done independently by 
Walker [64] and it is included in the COPT program. For the discrete case, a 
FORTRAN program is developed to carry out the estimator gain iteration until 
it approaches the steady state value. This program is validated by comparing 
the results of some examples with uncorrelated disturbance and noise run on the 
program DISC [28]. 

The experimental system has an unobservable sensor location for the 
accelerometer. It is verified to be the center of percussion (CP) [29]. This same 
physical behavior is used independently by Chiang in his mini-manipulator design 
[30]. Several different investigations of the experimental system observability by 
examining the observability matrix, the transfer function pole-zero loci and the 
modal disturbance and measurement distribution matrices are performed and 
some conclusions are given. 
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Finally, a gradient search scheme is developed to find the optimal sensor 
location based on the sensor location index which is similar to the performance 
index for the estimator design. This is a general scheme which is capable of solv- 
ing multisensor problems. A computer program is developed to do the search 
and a result for a single sensor case is given. The sensor location index (SLI) is 
then compared with the degree of observability (DO) suggested by Vander 
Velde. A double sensor example is given to conclude this chapter. 

Chapter III describes our experimental apparatus. There is not too much 
difference from that used by Rock {25] or Stoltz [26], except a new and thicker 
wing model is built to accommodate those accelerometers used as measuring 
instruments. A digital computer is used in this study for data collection and 
analysis. Some software is developed to facilitate the laboratory operation. 

Chapter IV describes the experimental methods and results. The approxi- 
mation to Timman's modified Theodorsen function suggested by Rock is reexam- 
ined. With the external sensors used and the doublet excitation he had found it 
impossible to verify experimentally this part of the model. However, it is a well- 
posed parameter identification problem [Refs. 31-32]. An on-line parameter 
identification scheme through Kalman filtering proposed by Mishne [33] is 
modified and used for this purpose. Simulations are carried out and found to be 
successful. But actual tests fail because of the existence of the process distur- 
bance and the measurement noise. 

A more accurate 5th order system model which includes those experimen- 
tally measured damping terms in the F matrix and takes into account the 
unbalanced force due to the applied torque in the G matrix. The dynamic cou- 
pling [Refs. 34-35] between the plunge and pitch modes through the actual 
wing section mass center offset from the elastic axis is used in the F matrix 
rather than using the distance between the entire plunge suspension system mass 
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center and the elastic axis. The accuracy of the system model is checked by exa- 
mining the theoretical and the experimental root locus versus varying wind 
speeds. The worst case for stable behaviors is about 1 % error in the pitch 
mode frequency. This same accuracy is also observed by comparing the computer 
simulated system responses to the linear motor doublet input with those of the 
actual system responses. 

The quality of estimation varies with the location and number of sensors. It 
can be indicated by the sensor location index we proposed. The difficulty in 
evaluating it as a function of the sensor location is presented in the estimator 
design section along with the approaches we use to get around it. The experi- 
mental test results are given and they prove the trends predicted by the theory. 

Chapter V gives a summary and recommendations for the future research. 

Appendix A includes the derivation of the system modeling and the 
unsteady aerodynamics model derived by Rock. Instrument calibrations and sys- 
tem parameter measurements are detailed in Appendices B and C. 

C. SUMMARY OF CONTRIBUTIONS 

1) The computer program OPTSYS is modified to the LOPTSYS to 
accept a more general input form including a rate measurement. The 
correlation between system disturbance and measurement noise is 
taken into account. 

2) An iterative scheme for calculating the discrete Kalman filter gain is 
derived to deal with the correlation between system disturbance and 
measurement noise for discrete systems. 

3) Experimental system observability is examined from several different 
approaches and some conclusions are given. System unobservable 
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sensor location is verified analytically to be the center of percussion of 
the two DOF typical section wing at zero wind speed. 

4) A tractable sensor location index is defined so that the optimization 
technique can be applied. A first order gradient search algorithm is 
developed and implemented. It is a general scheme which can be used 
for solving multisensor cases. The program doing the optimal sensor 
location search can report the figure of merit of the sensor location 
index for different sensor locations and increasing sensor number. 

5) An experimental apparatus is set up for feasibility studies. It includes 
a new and thicker airfoil with six accelerometers installed inside along 
the wing chord as measuring instruments. 

6) An on-line parameter identification through Kalman filtering is 
modified and used for aerodynamic parameter identification. 
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Chapter II 


PROBLEM FORMULATION 


A. INTRODUCTION 

State variable feedback control requires all the system states to be measured. 
However, the more frequent situation is that only certain linear combinations of 
the state can be obtained. 

By using a linear model of the system and the statistical models which 
characterize the system disturbance and measurement noise, Kalman developed a 
filter with which to reconstruct the missing states from any set of noisy measure- 
ments [Refs. 36-37]. 

Although a Kalman filter describes how to process the measurement data to 
achieve the optimal system state reconstruction, it doesn’t [38, p. 3]: 

1) solve the problem of designing in the presence of parameter 
uncertainties, 

2) provide a method for dealing with computational errors, 

3) solve the problem of establishing an optimal measurement schedule. 

The first problem is due to the presumption of a perfect mathematical model 
which requires the system be known exactly in order to formulate the optimal 
Kalman filter. Since this is often unrealistic and may lead to filter divergence, 
many authors have derived results to overcome this problem [Refs. 39-42]. The 
second difficulty is caused by the quantization limitation inherent in the com- 
puter implementation, especially in dealing with ill-conditioned problems. 
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Alternate recursive relationships have been developed, such as discrete square 
root filtering (Refs. 43-44], to cope with this difficulty. The third question is 
because there are various choices of the sensor type, number and location. 

The sensor type is usually limited by its utility, cost, availability, reliability, 
and other factors [45]; the sensor number is a tradeoff between the estimation 
quality and the cost of installing and maintaining, the weight as well as power 
consumption for extra sensors. All of them are highly cost dependent. However, 
a good choice of sensor location may increase the overall performance without a 
change in cost. Of course, the sensor location may still be constrained by its size 
and ease of access. 

While the sensor type is usually determined in the early design stages, the 
sensor placement is a degree of design freedom frequently not fully exercised, 
especially in large distributed systems. A standardized criterion could be very 
helpful to a designer in choosing the sensor number and their locations optimally. 
Some authors have already suggested approaches [Refs. 8-11], and the most 
intuitive result is given by Vander Velde [Refs. 12-14]. In this chapter, a new 
approach based on the steady state solution of the optimal estimation problem is 
presented to study the sensor location problem. The final goal is to develop a 
process that can be used to look for- the best sensor location directly and report 
its figure of merit for increasing sensor number to help the designer’s judgement 
in making his decisions. 

A experimental study for choosing the sensor number and their locations is 
conducted on a simple two DOF aeroelastic system. A detailed description of 
this experimental apparatus is given in Chapter III. Basically, it is a physical 
realization of an ideal two DOF typical section wing [Refs. 17-18]. The model- 
ing of the unsteady aerodynamics of this system was already developed and 
experimentally verified by Rock in 1978 [25]. So there is obviously an 
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advantage in using a similar system to save some time and effort. Besides, this is 
a simple physical system which is still capable of describing such complex aeroe- 
lastic phenomena as divergence and bending-torsion flutter thus making it a real- 
istic application area of our research objective. Also, many experimental studies 
have used this typical section wing, it facilitates the comparison of new and old 
experimental results [25, p. 2]. 

In this chapter, an example is first given to illustrate the effects of sensor 
types and locations upon the system observability. Secondly, the Kalman filter 
gains for both the continuous and discrete cases will be reformulated under the 
condition of correlated process disturbance and measurement noise in order to 
minimize the performance index chosen to be the weighted sum of the mean 
square estimation errors. Thirdly, a study of the observability of our experimen- 
tal system is investigated through several different approaches. The purpose is to 
find a direct connection that can be used to help determine the sensor location 
and number. Then a gradient search scheme is developed to achieve this goal. A 
comparison between Vander Velde’s criterion and ours is conducted assuming 
there is no process disturbance and using finite estimation time. Finally, a dou- 
ble sensor example is given to conclude this chapter. 

B. OBSERVABILITY OF DYNAMICAL SYSTEMS 

To many readers, the concept of observability of dynamical systems may 
seem to be a weak abstract idea rather than one with strong physical feeling. A 
real example showing its physical meaning should be useful for further investiga- 
tions. 
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Let us consider a uniform, free-free beam undergoing a transverse vibration 
as shown in Fig. 11*1 without taking into account the effect of gravity. Neglect- 
ing the shear distortion and rotary inertia, the equation of motion of the beam is 
governed by the Bernoulli-Euler model (46, p. 135]: 


d i y{x,t) a d 2 y(x,t) _ 

dx* El dfi 


( 2 . 1 ) 


where y[x,t) is the transverse displacement at station x and time t, and El, 
a , and L are the beam bending stiffness, mass density per unit length, and 
beam total length respectively. 



Fig. n-1 Uniform, Free-Free Beam Undergoing A Transverse Vibration 
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The natural boundary conditions are [46, pp. 163-166] 




(2.2a) 


>■« - - " 


(2.2b) 


Since the meaning of observability is understood most clearly in the modal 
coordinate system, the displacement is thus expressed as 


00 

fa*) = E (2-3) 

i=0 


where <f> t {x) are the mode shapes and q t {t) are the modal amplitudes [1, p. 36]. 

This separation of variables leads to a solution by substituting Eq. (2.3) 
into Eq. (2.1), and we obtain 


dx 4 




(2.4) 


and 


EL _L 

° <t>i dx 4 


_ 1 _ £0i_ 

<7,- dP 


= w?. 


(2.5) 


So the differential equations to be satisfied become 


dx 4 



(2.6) 
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and 


$ + ^ 


where 



After enforcing the boundary conditions, Eq. (2.2), we obtain 


(2.7) 


( 2 . 8 ) 


M x ) — A), A — ®» 

0i(x) = A t (x - ~), /?, = 0, 

and 

<f>,{x) = A{ [ (cos fcL - cosh/?,L)(sin/?,x + sinh/?,*) 
- (sin#,£ - sinh/5,L)(cos^,x + cosh/?,x) ), 

where /? satisfies 

cos/? t Lcosh^,L =1, « = 2, 3, • • • . 


(2.9a) 

(2.9b) 


(2.«c) 


( 2 . 10 ) 


The first four modes are plotted in Fig. H-2. The two rigid body modes 
correspond to the rigid body translation and rotation respectively. Since there is 
no node in the rigid body translation mode, a displacement sensor can always 
observe this mode. However, if we place a displacement sensor at those nodes of 
the other modes, i.e. at points B, D, E, G, I, the corresponding mode will 
not be observed at those particular sensor locations. This can also be seen from 
the exact pole-zero cancellations in their transfer functions (see discussions in 
Section D-3). If an orientation sensor is used and placed at those locations where 
the mode shapes have zero slopes, i.e. at points C, F, and H, 


- 14 - 



the corresponding mode will be missing from its measurements, that is, the orien- 
tation sensor can not detect the mode at a point of zero slope [47, p. 6]. This 
should provide a clearer picture of the dependence of observability on sensor 
types and locations. 



Fig. II-2 First Four Modes of the Transverse Vibration 
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C. KALMAN FILTER GAIN 


The system equations of motion of our experimental apparatus are derived 
in Appendix A, see Fig. A-l. They are presented here for ease of reference. 

• Equations of Motion without considering Aerodynamics: 



■ 0 

1 

0 

« 

0 


K 



2 w* 


2fa w or^2*a 


K 

K 


K 

K 

K 

K 


K 

a 


0 

0 

0 

I 


a 

• • 

a 



2 <*>* *a 

t j 2 

2? a W a 


a 



7* K 

ft 

7 2 K 

~ 

K 

4 


L 


0 0 

Kf l K 2 7 a b k t 

MfKb ' M T Tb I ar Kb 

0 0 

Kptg ._ 1 _ * a » K T 

M T 7*Kb 'la M t T 7%b 2 ' K 


Fv 

T v 


( 2 . 11 ) 
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• Equations of Motion considering Aerodynamics: 



( 2 . 12 ) 

where 

h, h, h: are the normalized plunge displacement and its time 

derivatives, 

a, a, a: are the pitch displacement and its time derivatives, 

Fy, Ty>. are the input voltages to the linear and angular 
motors respectively. 
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All the coefficients and constants are defined in Appendices A and C. 


C-l Continuous System 

Most of the mathematical representation of dynamic systems including the 
ones above can be described by the following linear, time-invariant, state 
differential equations: 


i(t) = Fi(t) + Gv(t) + r«<<), 
2 (f) = H4t) + JX*) + t<<), 

where 


z(<): 

«,X1 

state vector, with 


E { 4'o) 1 = E { M<oHb]M<oMW] T > = 

«(/): 

« t Xl 

control vector, 

2(f): 

n oftX 1 

observation vector, 

u<<): 

TlgXl 

white process disturbance vector, with 


E { u(f) } = 0, E { uMw t {t) } — QJ (f-r), 
or it can be simply written as v^t): N ( 0, Q a ], 


where Q v is the power spectral density of v\t), 
o(f): n ot Xl white measurement noise vector, with 

£{«(()}= 0, E { vft)v T (r) } = R,6 ((- r), 
or it can be simply written as «(f): N [ 0, R v J, 
where R v is the power spectral density of v(t), 
t t > <o time. 


(2.13a) 

(2.13b) 
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Furthermore, it is assumed that u {t) and v(t) are correlated, i.e. 
E { w(t)v T (r) } = V6 ( t-r ), where V is the cross spectral density of w(t) and 
v(t) [48], while the initial state z^) is uncorrelated with v{t) and v(/); u(<) is 
a given known input to the system. 

Suppose that a full-order estimator of the form [15, p. 351] 

40 = Fi{t) + Gu(0 + K(t)[ 2 {t) - H%t) - J u v(t) ] (2.14) 

is used to estimate the system state, and the estimation error is given by 

40 = 40 - 40- (2.15) 

The estimator can be designed optimally by finding the estimator gain matrix 
K(t), t 0 < r < t, and the initial condition iffo), so as to minimize the 
weighted sum of the mean square estimation errors 

J° = min E { e T (t)We(t) } = min tr [ WP{t) ], (2.16) 

m m 

where W is a symmetric positive-definite weighting matrix, and 
P(t) = E { e(t)e T {t) } is the estimation error covariance matrix. 

First subtracting Eq. (2.13a) from Eq. (2.14) and using Eq. (2.13b), we 
obtain 


e(t) = [F-tf(0tf]40 + tf(040-ri40, (2.17a) 

4*o) = e o = 4*o) - 4*o). (2.17b) 

It follows from the above equation that P(t) satisfies the following 
differential equation 
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(2.18a) 


m = [F-mum+miF-mn] 7 

+ Ki^Rjc^t) + rQ tt r 7 '-rvK r (o - K(t)v T r T , 

P[b) = Po- (2.18b) 


Substituting Eq. (2.18) into Eq. (2.16), and using the gradient table [49] in 
Appendix D to take derivatives, we obtain 


J = tr 


W 


/{ [F-K(r)H\P[r) + P[r)[F-mH\ T 

to 

+ K(t)R,K T (t) + TQJ t 
-TVK t (t)-K(t)V t T 7 *} rfr+P 0 j|, 


dJ 

dK 


= W I / { - 2P{t)H t + 2 K{t)R v - 2TV) dr 


f - = 2Wt R,dT > 0. 

9K 2 t 


(2.19) 

( 2 . 20 ) 
( 2 . 21 ) 


In order to minimize J with respect to K, 


dJ 

dK 


must vanish, that is, 


- 2 P{t)H t + 2 K(t)R v -2TV = 0, f>r><o, (2.22) 

which leads to 

K°(t) = [/W^+rKJf?; 1 , t>t 0 . (2.23) 

It can be seen that K°{t) b independent of the choice of the weighting 
matrix W. 
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Using the optimal gain matrix of Eq. (2.23) in Eq. (2.18), P(t) now 
becomes the solution of the matrix Riccati equation 

P(t) = [F-rVR-, l H]P{t)+ P{t)[F-rVR- v l H] T 
-P{t)H T R; l HP[t) + rQJ' T 
-TVR; l V T r T , t>t o, 

m = p» 

Since Eq. (2.16) can be rewritten as (15, p. 342] 

E{ e T {t)We{t) > = e T (t)We[t) 

+ E{[ e(t) - e(t) ] T W[ e(t) - Z(f) ] }, (2.25a) 

where 

e(t) = E{ e(t) }. (2.25b) 

It is a sum of two positive quadratic terms, the first term of which is obviously 
minimal when e(f) = 0. This can be achieved by letting e(< 0 ) = 0, because 
e(<) obeys the homogeneous differential equation 

2(/) = [ F- K{t)H]I(t) t (2.26) 

so the initial condition of the estimator can be chosen as z(<q) = % 

Since the system is time invariant, i.e. F, G, T, H, J u , Q w , R v , V are 
all constant matrices, it can be proved that the estimating process will reach a 
steady state as t — ► oo [16, p. 366]. 

The OPTSYS computer program developed by Hall and Bryson [27] uses 
the eigenvector decomposition technique to solve for the steady state estimator 
gain and root mean square (rms) estimation error. Since the use of accelerome- 


(2.24a) 

(2.24b) 
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ters as measuring devices is common in most control applications, OPTSYS has 
been modified by the author to the LOPTSYS to accept the more general input 
form 


i(i) = Fx(t) + Gu(t) + Tw{t), (2.27a) 

4t) = H zd i{t) + H^t) -f J' n u(0 + f v tt<f) + v' (t), (2.27b) 

where 

«<<)• N(0 ,Q a ], v'(t): N [ 0, R v ' ], 

E{ v{t)v' t (t) } = 0. (2.27c) 


Substituting Eq. (2.27a) into Eq. (2.27b), it can be shown that [Refs. 50, 
pp. 2-3, and 51, p. 73, and 52] 


00 = [ H :i F + H,\4t) + \ + K HO 

+ |H,jT + 4 M0 + »'(0, (2-28) 


so comparing with Eq. (2.13b), we have 


H = H zd F + H v 
J n = H zd G+f u , 

= H Z JT + J' v , 

tit) = lff zd T+ J' v Ml) + v' (t), 

r v = [//,/ + 4 m h,f + j’ w \ T + R v ' , 
v = Q W \HJ + J' W ) T . 


(2.29) 


It is important to know that the derivative term introduces a direct feed 
through of the process disturbance into the measurement noise and it is highly 
sensor location dependent. 
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C-2 Discrete System 

In this research the evaluation of estimation quality for different sensor 
numbers and locations is carried out on a digital computer. So the previous 
development has to be extended to the discrete case. The discrete estimator gain 
found in this section will be used for our discrete estimator implementation. 

First Eq. (2.28) is rewritten as 

40 = Hj{t)+ J u v(0 + J^t) + v’ (#). (2.30) 

Instead of modeling u(0 and v' (t) as white noises, we must think of them as 
colored noises with variances Q w and R n > , and correlation times T w and T v ' , 
respectively. That is, 

T 9 MO + MO = VM (231a) 

7>' (<) + *'(/) = ft AO, (2.31b) 

where rj w (t) and tfAO are white noise processes with power spectral densities 
2 T W Q W and 2 T v > R v ' , respectively. 

Let x a = £ x T w T v' r j and ij = £ »/ £ J , Eqs. (2.27a), (2.30) and 
(2.31) become 


*a{0 = F a x a(0 + Gv{t) + T^t), (2.32a) 

40 = H a x a (t) + JAO, (2.32b) 
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where 


[f r o 



H, = [tf J, /], 

9 = N[0, Q\, Q = diag[2r.9„ 2r, F?, ], 

which is a system with no white noise measurement. 

Thus the discrete system is of the form [Refs. 53-55] 


(2.33a) 


(2.33b) 


(2.33c) 


(2.33d) 


*«(*+!) = F d x ai k ) + GX*) + 

z(k) = H a x a (k) + JAk). 


(2.34a) 

(2.34b) 
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If we use a zero-order hold for the control, that is, 


«(') = <‘k), lk<t< 'hi. *= 0, 1, 2 (2.35) 


and assume that the observation is made at the controlling instant t*, then 


4+1 

*«(<*+i) = *(h+vh)*Jih) + I + l(h), 

4 


Ah) = H a x a {t k ) + JXh)- 


where $(f*+i,f*) is the state transition matrix from time t k to t k+l . 
Comparing Eq. (2.34) with Eq. (2.36) we obtain: 


Fi = $(**+!>**)> 

<i+i 

Gd = / ${h+v T )GdT, 

4 

4+i 

VAQ = V(h) = / *(**+i> r )r .«?(»)*, (2.37a) 

4 

and 

E {*„(*«) > = E{ W*b)- ? JW*h)- ? J r > = Mo, 

£ { ui*) } = 0, <?, = £ { vMnl(t ) } 

<i+i 

= I m + iJ)r a Qr^ T (t k+1 ,r)dT. (2.37b) 

4 


In the usual case in which the sampling instants are equally spaced, i.e. 
t k+1 - t k =' T t , while the system is time-invariant, then 


(2.36a) 

(2.36b) 
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Ft = *(w») = = e F - T -, 

T, 

G t = ( fe F -’dr\G, 

0 


T, 

Ql = / /•' r.QrJe^Vr 




(2.38) 


Similar to the continuous case, we can design an optimal discrete estimator 
as [15, p. 550] 


x a (k+l) = F a (*+1) + K(k+1) X 

[ z{k+ 1) - H a x a (k+1) - J u v(k+ 1) ], (2.39a) 

where 

*„(*+!) = F d xJ[k) + G d u(k), (2.39b) 

by finding the gain matrix K(j), kq < j < k, and the initial state i a (A^) in 
order to minimize the weighted sum of the mean square estimation errors 

J° = min E[ e T {k)We(k) } = mm tr ( 1 VP[k) ], 

K(k) K[k) 

where W is a symmetric positive-definite weighting matrix, and 
4k) = i,(k) - xjk), P[k) = E { e(k)e T (k) }. 

Subtracting Eq. (2.34a) from Eq. (2.39a) and using Eqs. (2.34b) and 
(2.39b), we obtain 


(2.40a) 


(2.40b) 
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4*+l) = I /- 1 1 7„(*+l) - xtk+1) ]. 


(2.41) 


It follows from the above equation that P{k) satisfies the following recurrence 
relation, 

P[k+l) = [ I- K(k+l)H a ]M{k+l)[ I- K(k+l)H a ] T , 

where 

m+l) = E { ( iji+l) - xjk+l) ][ 7,(*+l) - xjk+l) ] r } 

= E { [ Fje(t) - ri/k) ]( F d e{k) - y/k) ] r } 

= FAk)Fl + Q i} 

which is the covariance matrix of the estimation error before measurement. 
Substituting Eq. (2.42a) into Eq. (2.40a), J can be found as 

J = t r { w[[I-Wk+l)H a ]Wk+l)[I-K(k+l)H a ]^}, (2.43) 


(2.42a) 


(2.42b) 


and its derivatives are 

W [ - 2( /- K(k+l)H a )M(k+l)H ^ , (2.44) 

2W[ H a M(k+l)Hl\ > 0. (2.45) 

dJ ' 

J can be minimized with respect to K by letting -^77 equal to zero, that is, 

dK 

- 2[ I- K{k+l)H a ]M{k+l)Hl = 0. (2.46) 


dJ 

dK 

d 2 J 
dK 2 
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From the above equation, K°(k+l) can be solved as 

K'(k+1) = (2.47s) 

with initial condition 

Ai(*b) = Mq. (2.47b) 

Substituting Eq. (2.47a) into Eq. (2.42a), then 

tf(*+l) = ( /- K°(k+l)H a ]M(k+l). (2.48) 

As in the continuous case, we can let e(J^) = 0 by choosing 

= W (2.49) 

For a time invariant system, i.e. F * G dt H a , J w Q d are all constant 
matrices, it can also be proved that the estimating process will reach a steady 
state as in the continuous case. A computer program DISC developed by Katz 
[28] can be used to find the steady state solution, but it is not able to handle the 
correlation between the process disturbance and measurement noise. So a pro- 
gram coded in FORTRAN was developed by the author to iterate Eqs. (2.47- 
2.48, 2.42b) until the steady state in Eqs. (2.47-2.48) is reached. Better tech- 
niques can be used to find the steady state solution at faster convergence speed 
but they are not considered here. 
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D. EXPERIMENTAL SYSTEM OBSERVABILITY ANALYSIS 

For ease of analysis and notation, we will use the continuous formulation 
hereafter, and the system is assumed to be time-invariant. 

D-l Differential Equation Investigation 

As derived in Appendix B-l, the system equations of the physical system, 
Eqs. (A-23)-(A-24), are repeated here 

T 

MjH + S Q d + C k h + K h h = F j- - L, (2.50a) 

S Q h + I Q d + C a a + K a a = T a + M a . (2.50b) 


Without considering the driving forces and aerodynamic forces, i.e. in vacuo, 
and neglecting the damping effects, the equations normalized by the the semi- 
chord can be written as 


r 4 r , K 2 u 2 a z a 

* = 'T h+ ~~K~ a ’ 


(2.51a) 


a = 


r <4 


K? a 




(2.51b) 


where all the coefficients are defined in Appendix C. 

Using a linear accelerometer located at a distance d away from the elastic 
axis (EA) for linear acceleration measurement, 7, see Fig. II-3, 

7 = h + da 

_ oj\ dx _ 

= —■( K 2 x q -d)a + -£(-^-l)h, (2.52) 

•Vv 
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and the second derivative of this measurement is 


, <4<4*« , ^ _ Tx «1 , . % , r- 

2 1 ( A ’ A* ( ;* 11 

<4 _ Knu&JLx n dx~ 

+ l~^(K2*a-d) + - J ^( 1 ?-- Dior. 


(2.53) 


Combining Eqs. (2.52) and (2.53), we can solve for h and a in terms of Jand 
z provided the determinant 


A 


44 

K? 


[(K 2 x a - J)u£-( 



44 . 

K 2 1 



? + ( <4 - < 4 ) <* - *& 4 1 


(2.54a) 

(2.54b) 


does not vanish. 



Fig. II-3 Acceleration Measurement 
(An accelerometer located at x, w.r.t. midchord (MC) 
is d from the elastic axis (EA)) 
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Determinant 


Since d — x t - a, the determinant A can be plotted as a function of x t 
for given system parameter values: 

— 1821.58 (rad /sec ) 2 

u> 2 = 3737.83 (rad /sec ) 2 

K 2 = 0.3684 

AT = 0.96469 

x Q = 0.19 

d = x t - a 

f a =* 0.37668 

a = -0.3 

: - 1. < ?, < 1. (2.55) 


as in Fig. 11-4. 



-10 - 0.6 - 0.2 0.2 0.6 1.0 


— Sensor Location x s 

Fig. II- 4 Determinant A va Different Single Sensor Locations x, 

It can be seen that A = 0 at x t - 0.2, that is, the system will lose its 
observability if we put an accelerometer at this location. The following example 
will be used to explain the physical meaning of this singularity. 


- 31 - 




Example: (see Figure II-5) 

Suppose an airfoil with mass, m and moment of inertia about EA, 4 ■ 
suspended at point O, which is of distance d away from EA. Assume a linear 
spring with spring constant K k and an angular spring with spring constant K a 
is attached to the airfoil at EA. Now displace EA laterally away from its neu- 
tral position a small distance x and release, the reaction force at point O in i 
direction R g can be solved as follows. 



Fig. n-5 Example 

The force and moment equations about point O are 


R s -K h x = mx e = m(d-x a )0, 


K a - d - K k xd = I e 6, 


so 


B, => — 7 -- * [ -K 0 j - K k ,i) + K t i. 


- 32 - 



Since 


K k = tnul 
K a = I a ul 

I 0 = mi 2 = mi 2 - mfi+nid- x a ) 2 , 

and 

= 4-4 + ( rf - x a) 2 — r l~ 2 x a d + <**, 
so 

Rj = --- - [ -mJlxd- mi 2 u%-2- ] + mw|i 

K d 

= ^ " 2* a d + d 2 ) - (d - z a )dwj[ - 1 2 <J* J d ^ ] 

= 4 1 ^ ^ )• < 2 - s6 > 

*0 *ci 

Thus if d is so chosen that R x = 0, the accelerometer will have zero 
measurement at that position. Physically, this is the center of percussion (CP), 
where there is zero specific force. Comparing with Eq. (2.54a), they are 
equivalent when K 2 — 1. Since plane motion is involved in this example and 
when there are no other forces applied, the instantaneous center will be located 
at the center of percussion, where the ratio of linear and angular accelerations is 
such that there are no movements during vibration [29, p. 481). When external 
forces are applied, the acceleration measurements at CP may not be zero but 
the ratio of the incoming linear and angular components is such that there will be 
no distinctions between them. This can be seen by letting 

< 4 (*. - <0 = - 1 ) 

n 

in Eq. (2.52). Then the measurement contains only one of the frequency com- 
ponents, either plunge or pitch but not both. 
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D-2 Observability Matrix 

A more systematic way to check the observability is to write the system 
equation in the state space form as Eq. (2.13). Equation (2.11) is first used for 
this investigation. If only one accelerometer is used, the measurement equation 
can be written as 

it) = HJW + tit), 

where 

— [° 1 0 (*,-«)], 

and 

it) = (A(<) ^*) a(t) a(t) ] T . 

Using Eq. (2.29), we have for our physical system that 

H = H zi F, 

Ju = H zd G, 

it) = HJxit) + v ( t ), 

V = QJ T Hl, (2.58) 

where we assume V =0. 

Then the observability matrix O gf) can be constructed as 

H 
HF 
HP 
HP 




(2.57a) 

(2.57b) 

(2.57c) 
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In general, O ob is an (n o6 X n,) by n, rectangular matrix. The deter- 
minant of the n, by n, square matrix 0j b 0 ob , det (Oj b O ob ), can be plotted as 
a function of x t . In this formulation the aerodynamic effect can be considered 
easily by using Eq. (2.12) which adds an additional aerodynamic state. Now the 
measurement distribution matrix becomes 


Hu = I o i o (?.-«) o i, 

4‘) = [*<) R>) a(Q MO H01 T < 


and the observability matrix will be 



- 

H 

HF 

HF* . 

HF* 

HF* 


(2.60a) 

(2.60b) 


(2.61) 


The det (Oj ft O ofc ) can thus be plotted for various wind speeds. Figure H-6a 
is a plot of the theoretical estimation error performance index tr ( WP) versus 
different single sensor locations x t at zero wind speed, where 
W = diag [ 1, 10 _1 , 10, 10" 2 , 1 J. This result is obtained from the steady state 
covariance analysis, it gives finite estimation error for finite power spectral den- 
sity inputs of the process disturbance and measurement noise even at x #m i n . 
Figure II-6b shows the det ( Oj b O ob ) versus different single sensor locations x t 
at zero wind speed. For different wind speeds, the plots look the same except the 
x $ min corresponding to the minimum value of det (0^ b 0 ob ) will be different. 
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Fig. II-8b Determinant det {oLO oi ) 
vs Different Slagle Sensor Locations x, at Zero Wind Speed 
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By comparing the plots in Fig. U-6, we can see that x t min corresponds to 
the unobservable location of the sensor. The same conclusion also holds for vary- 
ing wind speeds. Otherwise, the absolute magnitude of det (Oj b O ob ) bears no 
direct relationship with the observability of the system. This can be seen that 
for the same magnitude of det ( Oj b O ob ) the corresponding sensor locations can 
have different results for tr ( WP). Actually, the det ( Oj b O ob ) can only provide 
binary information about the system to be analyzed, i.e. either observable (when 
det {Oj b O ob ) 0) or unobservable (when det {Oj b O ob ) = 0) [16, p. 457]. 

D-3 Transfer Function and Modal Form 

Using the state space formulation Eq. (2.13a) and the measurement equa- 
tion Eq. (2.13b), we can find the transfer function from u^t) to 40 38 


As) 

W{s) 


8Hi\ 8i-\)~ l r l r, 


(2.62) 


and the modal form as 


£(t) = A£(0 + r l Gv(t) + r l Tw{t), (2.63a) 

40 = HT{(0 + Jn<0 + <0, (2.63b) 

where A = T~ l FT, f(<) = T' l 40> aQ d T is the right eigenvector matrix of 
the open-loop system. 

The product of disturbability times observability can be determined by exa- 
mining the transfer functions or the modal distribution matrices. Figure II-7 
shows the pole- zero loci for different single sensor locations and various wind 
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speeds. The conclusion from this figure is that at the unobservable or undisturb- 
able position of the sensor, there is an exact pole-zero cancellation of the pitch 
mode. It seems that as the separation between the pole and zero increases, the 
estimation error for that mode decreases. Although this qualitative behavior is 
general for all cases, we did not use it for best sensor location searches. 

If ( z , w) were scaled so that one unit of each component is of comparable 
significance, then {HT); and (7*‘ 1 r) l - indicate relative observability and distur- 
bability of «th mode. Outer product of these two vectors forms a rank one resi- 
due matrix. The norm of this matrix equals the product of the length of these 
two vectors, and it is a measure of the significance of »th mode in the input- 
output relation. 

A thorough examination of the residue matrices reveals useful information 
about its connection with the magnitude of the estimation error and it can help 
to locate the unobservable or undisturbable location of the sensor. 

E. CRITERIA FOR CHOOSING THE BEST SENSOR LOCATION 

The choice of sensor location represents a significant design freedom avail- 
able to the designer. However, it is usually not a very straightforward choice to 
make. To make a good choice, it is necessary to have a standardized criterion to 
choose the sensor number and their locations. Most important, it must be a 
quantitative measure that can be easily computed and have a physical interpreta- 
tion so that engineering judgement can be based upon it. In order to get a per- 
spective on this research, it may be helpful to review the criterion proposed by 
Vander Velde [12-14]. 
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E-l Degree of Observability (DO) 

For the Kalman filter with uncorrelated process disturbance and measure- 
ment noise, i.e. V = 0, from Eq. (2.24) the estimation error covariance matrix 
can be found satisfying the following equations 

P[t) = FP{t) + mF r -P{t)H T R; 1 HP{t) + rQJ' T \ <><o, 

(2.64a) 

P(to) = P 0 . (2.64b) 

Vander Velde chooses Q w = 0 by arguing that it is an external influence 
not related to the sensor set (however, this is not true when accelerometers are 
used as we pointed out at the end of Section C-l) and uses the information 
matrix S(t) — formulation, then Eq. (2.64) becomes 

S(<) = -S[t)F-F T S{t) + H T R; l H, t>t 0 , (2.65a) 

S(*o) = S 0 , (2.65b) 

using the fact that 

m = - wmm w 

Assuming there is no information about the initial state, i.e. P 0 = oo or 
S 0 = 0, one b interested in measuring how much information has been accumu- 
lated in S[ T t ) up to a specified estimation time T e , that is, the size of 5( T e ). 
This can be obtained by reference to the quadratic surface in the x-space, i.e. 

z T S~ i x = 1. (2.67) 
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After the transformation 


V = Tx, 


( 2 . 68 ) 


where T = diag | e lnMt . . . , | . . . , | | J and are the max- 

imum errors that can be tolerated in the directions x i} the volume of the 
transformed ellipsoid can now be measured in the equimeasure space. The degree 
of observability (DO) is defined as the maximum spherical volume (V$) con- 
tained in the ^-dimensional ellipsoid ( V E ) plus a lesser weighted excess volume 
( Vg- V s ) due to the nonideal volume distribution. Specifically, 


DO 




with 



(2.69a) 



(2.69b) 


and Vi are the eigenvalues of TS[ T e ) T. 
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E-2 Sensor Location Index (SLI) 

In Section A, the performance index for optimally designed estimator at 
steady state is found to be 

J° t , = min tr [ WP„ ], (2.70) 

Am 

where P ge is the solution of the steady state Riccati equation 
[ F-YVBr,'H\P„ + P„{ F- TVR?H\ t - 

+ yqj t - rvR-„ l v' r r r = o. (2.71) 

It can be seen that P sl as well as J° tt are functions of the measurement distri- 
bution matrix H. Specifically, they are functions of the sensor location along the 
wing chord, x t , for single sensor case. This is also true for multiple sensors, in 
which case, x t will be interpreted as a vector. The performance index can be 
plotted as a function of x„ and by a simple judgement, the location where it 
has the smallest value should correspond to the best sensor location under given 
conditions. 

The above conclusion suggests that a similar sensor location index 

J°sl = min J°„ = min { min tr ( WP tt ) } (2.72a) 

H H Aji 


can be used to find the best sensor location, where W is a symmetric positive- 
definite weighting matrix and can be chosen similarly to T as 
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To incorporate the constraint of Eq. (2.71), we introduce the use of a 
Lagrangian multiplier matrix A, and denote the Eq. (2.71) by Then 

Eq. (2.72) can be written as 

J°SL — min tr ( WP„ + A AP„,H ) ), (2.73) 

H 

since 

A AP">H) = 0- (2-74) 


After substituting and using the gradient table in Appendix D and Eq. (2.23), 
we can find that 


1^2- = [ F- K,J1\ T S. + A| F-K„H] + w, (2.75) 

dJqr v 

~ = - 2 KfrP... (2.76) 

The gradient of J SL w.r.t. each individual sensor location is obtained by apply- 
ing the chain rule. Using Eqs. (2.58) and (2.60a), we have in our case that 

= - 2 KfrP.fl, (2.77) 


where F 4 y is a 1 by n, matrix which is identical to the 4th row of the F 
matrix. For multiple sensor case, x t can be considered as a n ob by 1 vector. 

dJ SL 

In order to minimize Jgi, both — — and — — must vanish. 


dP, 




dH 
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E-3 Gradient Search Procedure 

The gradient searching process for the optimal sensor location can be pro- 
ceeded as: 

1. guessing initial value of x e with which to compute the component 
values of the H matrix, 

2. using LOPTSYS to solve for P tl and K t „ 


solving the Riccati equation of Eq. (2.75) for A, 

dJct <9./c/ 

computing the gradient — of Eq. (2.76), and — — of Eq. (2.77), 

OH dx t 

9J S l 

if — — < € (a small number close to zero given by the designer), 


stop; otherwise 


6. computing 




dJ SL 

dx, 


A x, 

•i 


(2.78) 


under the constraint of -1 < x, < 1, 

7. go back to step 2. 

This procedure is also shown by the flow chart in Fig. II-8. 

Figure II-9 shows the sensor location index J SL versus different single sen- 

_ dJsL 

sot locations x„ and Table II-l gives the values of J SL and — — . From 

dx e 

dJ SL 

which we can see that — 3- truly reflects the slope of the curve and a local 

dx, 

minimum exists. 
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Fig. II- 8 Gradient Search Procedure for the Optimal Sensor Location 
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Using the gradient 


for sensor location search is easy and convenient, 


dJsL 

dH 


however, 


dJ< 


SL 


dH 


= 0 


can lead to a local minimum. 


So the use of it must be 


with care. 


E-4 Comparison between the SLI and DO When Q v — 0 

Since in the actual applications, the process disturbance may be close to 
zero, i.e. Q w — *• 0, and V -* 0, then the P tl will approach zero. In this case, 
we could take into account the estimation time T t to avoid the singular result 
and solve Eq. (2.24) by integration. 

Figure 11-10 shows the result by letting Q w = 0, V = 0 and T e = 1 sec. It 
can be seen that the trailing edge is still the best location to install an accelerom- 
eter while the center of percussion is still the worst. 

Using Vander Velde’s information matrix formulation, Fig. 11-11 shows the 
sum of its weighted diagonal terms versus different single sensor locations. It 
shows the trailing edge is the best sensor location, however at the center of per- 
cussion the index value does not go to zero. Because, in modal form, only one 
diagonal term vanishes when the sensor is at the center of percussion. 

Figure 11-12 is a plot of the information matrix determinant versus different 
single sensor locations and it shows zero value at the center of percussion while 
the best sensor location is still at the trailing edge. 

For finite time estimation, the SLI is independent of the initial value 
assumptions because of their fast decaying rates for stable systems. However, 
this is not true for the DO. For some cases, the available initial information 
may dominate the value of the DO and makes the sensor location preference 
indistinguishable. 
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PERFORMANCE INDEX, trace(WP) 


WIND SPEED * 1.00 m/sec 



— t.O'Q: -0.75 -0.50 -025 0.00 0.25 0.50 0.75 1.00 


LOCATIONS OF SENSOR' 

Fig. n-10 Sensor Location Index J s r 
vs Different Single Sensor Locations x t for Q w = 0, V = 0, T e — 1 sec 
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PERFORMANCE INDEX, trace(S/W) 


WIND SPEED = 1.00 m/sec, ESTIMATION TIME = 1.00 sec. 



— 1.00 - 0.73 - 0.50 - 0.23 0.00 0.23 0.30 0.73 1.00 


LOCATIONS OF SENSOR 


Fig. n-11 Performance Index tr ( S/ \V) 
vs Different Single Sensor Locations x, for T e = 1 sec 
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PERFORMANCE INDEX, det(S) 


WIND SPEED- 1.00 m/sec, ESTIMATION TIME = 1.00 sec. 



LOCATIONS OF SENSOR 


Fig. 11-12 Performance Index det (5) 
vs Different Single Sensor Locations x, for T e — 1 sec 
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E-5 Relative Merits of the SLI versos DO 

Compared with the degree of observability (DO) suggested by Vander 
Velde, we have chosen to use the sensor location index (SLI) because it is more 
convenient for our purposes for the following reasons: 

1) Practically, it is much more meaningful to measure directly the state 
estimation error than the hypothetical ellipsoid volume, especially 
when the system state has physical meanings. 

2) The measurement noise will include the process disturbance when 
accelerometers are used, see Eq. (2.29). This inclusion of the process 
disturbance is highly sensor location dependent and it may affect the 
choice of the optimal sensor location. The sensor location index takes 
this into account in its formulation. However, in our case both the 
SLI and DO give the same answer. 

3) When evaluated for a finite time to avoid singular results the sensor 
location index is independent of the initial value. 

4) The existing design program LOPTSYS and some of its supporting 
subroutines were developed at Stanford. They are available and 
make it straightforward to compute J SL ) the value of J SL could be 
computed as a by product during the design process. 

5) Most important of all, this choice of sensor location index J SL makes 
the gradient formulation easy to handle and compute thus making 
J SL tractable for the optimal sensor location search. 
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F. DOUBLE SENSORS CASES 


Figure 11*13 shows the theoretical results of single sensor as well as double 
sensors, with one sensor fixed, versus varying sensor locations along the wing 
chord. Due to the existence of Q w and V, there is no simple relationship 
between the single sensor and collocated double sensors at the same sensor loca- 
tion. Adding a second sensor generally improves the SLI, with the improvement 
increasing as the second sensor moves toward the trailing edge. 

The best location for the second sensor is still at the trailing edge, same as 
the first one. This result shows however the negative side of our sensor location 
index definition, namely it is too dependent on the mathematical modeling of the 
physical system. For example, a practical estimator designer would rather choose 
to put one sensor at the leading edge and one at the trailing edge. Their meas- 
urements can then be used to solve for the linear and angular accelerations which 
are independent of the system modeling. So the two sensors can get 2 DOF 
information directly albeit the initial conditions. Since these dc information is 
rarely needed for flutter suppression, it suggests that a frequency-band weighted 
SLI may be useful to cope with system parameter uncertainties and modeling 
errors. A more realistic criterion which takes this into consideration will be 
recommended for the future research. 
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Chapter HI 


EXPERIMENTAL APPARATUS 


A. INTRODUCTION 

The idea of the typical section wing was devised during the 1930’s by 
Theodorsen and Garrick to simulate an actual wing by matching the properties 
at a station 70-75% of the distance from the centerline to the tip [18, p. 189]. 
Ideally, as shown in Fig. HI-1, it consists of a thin, rigid wing immersed in a two 
dimensional, incompressible airstream, suspended by two sets of springs allowing 
elastically restrained but uncoupled rotation about the elastic axis (EA, line of 
shear centers [18, p. 281]) and translation perpendicular to the free stream 
airflow. 

This chapter describes an experimental apparatus that is the physical reali- 
zation of the ideal two DOF (degree of freedom) typical section wing. Its 
suspension system was built by Rock in 1978 to investigate the aeroelastic sta- 
bility of a simple, two DOF wing in Stanford’s small, subsonic wind tunnel 
[25], This suspension system provides spring restraints for plunge and pitch 
movements while allowing virtually no elastic coupling between these two 
motions. It is also characterized by small structural damping and relatively high 
stiffness in the remaining DOF’s. 

A schematic view of the complete system is shown in Fig. IE-2. A descrip- 
tion of the major subsystems is presented, followed by a summary of system 
parameters and a discussion of the system performance. 
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Fig. m-1 Typical Section Wing 


B. WIND TUNNEL 

The wind tunnel used in this research at Stanford University is a small, 
closed-circuit, subsonic, vented tunnel, see Fig. IE-3 for its plane view [56]. It 
has a dynamic pressure range of 200 N/m 2 to 2000 N/m 2 (an airspeed range 
of 19 m/sec to 65 m/sec, M miX = 0.2). The dynamic pressure is controlled 
by adjusting the pitch angle of a 16-blade constant speed (19 Hz) fan. Five 
screens upstream reduce the mean turbulence level at the test section to approxi- 
mately 1 %. The contraction ratio is 8.67 and the total distance around the 
center line is 25.1 m. 


- 55 - 



TORQUE 

MOTOR 


Fig. m-2 Experimental Apparatus 
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Fig. HI-3 Plan View of the Stanford Subsonic Wind Tunnel 


The test section is 0.4572 m (18 in.) square by 0.9017 m (35.5 in.) long. It 
consists of a welded steel frame with three removable Lexan walls and one steel 
bottom wall, mounted on a welded steel cart with castors for ease of handling. 

The suspension system and actuators are bolted to mounts welded to the 
test section. All have been machined in place to provide proper alignment and 
rigidity. The only hardware inside the test section is the airfoil with its endplates 
as well as those sensors inside the airfoil for the “in-flight” estimation. All 
suspensions, actuators and monitoring sensors are located externally. 

C. AIRFOIL 

The airfoil is constructed of a fiberglass-laminate skin wrapped around a 
foam core. It has a NACA 0015 profile with a span of 0.4191 m (16.5 in.) and 
a chord of 0.2413 m (9.5 in.). A core of styrofoam is cut to the shape using a 
hot-wire passed over two aluminum templates (shown in Fig. HI-4) numerically 
machined to NACA 0015 specifications [57, p. 113]. 
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Four layers of bidirectional weave fiberglass cloth and epoxy cover the core 
to provide structural rigidity and one layer of light weight fiberglass cover pro- 
vides a smooth aerodynamic surface. Surface irregularities are removed by sand- 
ing and filling with quartz microspheres and resin. 

The wing spar is initially a single piece of square aluminum tubing, 19.05 
mm (0.75 in.) on a side and 3.175 mm (0.125 in.) thick, centered on the mean- 
line at 35% of the chord behind the airfoil leading edge and extending out each 
wing tip to the external suspension. It is cut into two pieces to fit in a rectangu- 
lar aluminum box which is constructed to mount an angular rate measuring dev- 
ice. The protruding ends on the box are fitted closely inside the cut spar tubing 
thus transmitting the structural loads across the wing, see Fig. HI-5. Six 
accelerometers of two different sizes and sensitivities are installed inside of the 
airfoil for the “in-flight” estimation. They are mainly laid out in the middle of 
the wing section to maintain the symmetry of mass distribution, and fixed onto 
the surrounding styrofoam and one side of the fiberglass skin through epoxy and 
microsphere filler. 

Circular endplates are fastened to each end of the airfoil, with a diameter a 
little bigger than the wing chord, machined from plexiglass 12.7 mm (0.5 in.) 
thick. A gap of 6.35 mm (0.25 in.) exists between each endplate and the wind 
tunnel wall to avoid boundary layer effects and maintain two-dimensional flow 
over the airfoil [26, p. 32]. Since both endplates move with the wing, they con- 
tribute significantly to the total apparent mass (10%) and moment of inertia 
(36%) of the system. 
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D. SUSPENSION SYSTEMS 


The suspension system is designed to constrain the airfoil to move in the 
plunge and pitch DOF’s without introducing any friction or elastic coupling into 
these motions. Hence, the suspension system employed can be treated statically 
as two independent systems — one for plunge and one for pitch. 

D-l Plunge Suspension 

The plunge suspension consists of four folded-cantilever springs mounted on 
the test section exterior, connected across the section top and bottom by light 
weight magnesium-aluminum cross beams. The cross beams are also connected 
along one vertical side of the test section by a third beam, as shown in Fig. III-2. 

This arrangement constrains the airfoil to the plunge motion only (Fig. III- 
6, direction 2), and it is stiff in all the other directions. The only suspension 
losses are in the material flexural hysteresis and are very small. The plunge 
direction is designed to be horizontal rather than vertical to avoid the gravity 
bias. 

The diagram of a single spring is given in Fig. HI-7. It is stiff in directions 
2, 3, and 4, but compliant in 5, 6, and 1, the desired direction. This spring can 
be analyzed as a group of four cantilever beams, two in parallel forming the 
center web, and two parallel split beams, each containing one half of both outside 
webs, which are in series with the central ones. The spring rate for displacement 
taken at the free end of a cantilever beam with the constraint condition of zero 
slope, is given by 


K, 


12 El 
/3 ’ 


(3.1) 
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Fig. Ill- 8 Arrangement of Plunge Suspension Springs and 
the Possible Degrees of Freedom of Motions 



Fig. m-7 Folded Cantilever Spring 


- 62 - 


from elementary beam theory [58, p. 212], where /= wfi/ 12. Since the four 
beams form a series-parallel set, the overall spring rate is the same as that of a 
single beam, given in Eq. (3.1). 

Each spring was machined from Beryllium-Copper (E = 1.28 *10 n N/m 2 ) 
with t — 1.68 mm, w — 28.6 mm, and l = 203 mm. The resulting spring 
rate is 

K , = 2075 N/m, (3.2) 

and the total spring rate along the plunge direction is 

K h = 4AT, = 8300 N/m. (3.3) 

The measured spring rate of the actual system is 10383 N/m, 25% higher 
than its predicted value due to the uncertainties in the parameters E, l , w, but 
principally in t. The springs are offset by 10 mm to avoid snap-through at 
their zero position. Adjustable stops on the top and bottom cross beams allow 
the plunge travel to be restricted if necessary. They proved to be useful when 
operating in dynamically unstable situations. With the pitch DOF locked (see 
D-2), the total effective plunge mass and the structural damping in the plunge 
mode were experimentally determined by using an FFT analyzer, see Appendix 
C for details. 

The vertical reinforcing beam shown in Fig. IE-2 helps stiffen the suspen- 
sion in the differential bending mode. The plunge actuator and plunge position 
sensor are located at the center of this beam. This location is chosen because it 
is at the node of the differential bending mode. Consequently, the actuator does 
not excite this mode and the sensor does not observe it either. 
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D-2 Pitch Suspension 

The pitch suspension system is formed by torsional springs fastened between 
housings on the wing spar tips and the plunge cross beams as shown in Fig. III-8 
for details. The springs are double cantilever Bendix flexural pivots, which allow 
elastically restrained pitch motion about the wing spar centerline. They are com- 
mercially available in a variety of spring rates and sizes. The bushings shown in 
Fig. DI-8 allow the housings to accommodate different sized pivots. The spring 
rates for the plunge and pitch suspensions are designed or selected to allow the 
test system to have an adequate separation between those two major system 
modes and the desired flutter frequency. The two flexural pivots used have a 
measured combined spring rate of 43.04 N-m/rad. Adjustable stops are also pro- 
vided for the pitch motion. 

Since the mass of the plunge suspension is elastically decoupled from the 
pitch rotation, the actual mass of the system rotating about the pitch axis is con- 
siderably less than the total system mass M T . The experimentally measured 
pitch mass Af 2 , including the wing section and its endplates and part of the 
torque motor, is 2.075 kg. Its center is 0.0165 m behind the elastic axis thus 
yielding a value of S Q = M 2 x a = 0.04814 kg-m. 

The effective moment of inertia of the wing section and the pitch mode 
structural damping were experimentally determined by using an FFT analyzer, 
see Appendix C for details. 
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E. ACTUATORS 


Two actuators are used, one for each degree of freedom. They may be used 
independently or simultaneously to simulate the effect of an aerodynamic control 
surface. 

E-l Plunge Actuator 

The plunge actuator is a linear motor (see Fig. HI-9) of the type found in 
computer disc drives. It has a bandwidth of 65 Hz and is capable of 62 N (11.7 
Vmput ) f° r short periods of time and 22 N (4.15 V input ) for continuous opera- 
tion. 

The voice coil current to force relationship is only linear within 5% but the 
repeatability ( 1 %) is such that the motor can be easily compensated. The 
motor is driven by a current-drive amplifier to eliminate the back electro- 
magnetic force (EMF) damping. 

E-2 Pitch Actuator 

The pitch actuator is an Aeroflex brushless torque motor shown in Fig. III- 
10. It is also driven by a current-drive amplifier. It has a bandwidth of 100 Hz 
and is capable of 2.12 N-m (9.5 Vin put ) for short periods of time and 1.06 N-m 
(4.75 V input ) for continuous operation. 

The motor (2.27 kg) acts through the four-bar linkage shown in Fig. HI-10 
to reduce its apparent mass added to the airfoil by roughly 100:1. 




Fig. HI-9 Linear Motor and LVDT with its Suspension 
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The operation of the linkage is explained with the aid of Fig. 111-11. The 
elastic axis of the airfoil is constrained to move rectilinearly by the plunge 
suspension, Fig. Ill- 11a. As the airfoil plunges, the linkage deflects without pro- 
ducing any torque (but it does produce an unbalanced force along the cross 
beam) about the elastic axis, Fig. Hi-lib. At any plunge position, however, the 
motor can transmit a torque through the linkage into the elastic axis, Fig. IH- 
11c. All joints in the linkage are very compliant flex pivots. There are no bear- 
ings used in order to reduce the friction effects and maintain the repeatability of 
the system performance. Consequently, the only suspension losses are in the 
material flexural hysteresis and can be neglected. 



Fig. m-11 Operations of Four-Bar Linkage 
a). Normal b). Translation c). Rotation 


- 88 - 


F. SENSORS 


A monitoring sensor is provided for each degree of freedom. Six accelerome- 
ters inside the airfoil can be used in various combinations for the “in-flight” 
estimation. A pitot-static tube is provided for airspeed measurement. 

F-l Plunge Sensor 

The plunge displacement transducer is a Schaevitz Engeering type 500 HR 
LVDT (linear variable differential transformer) which is excited with a 400 Hz 
signal, and when demodulated and filtered, produces an output voltage propor- 
tional to the plunge displacement. It is linear within 1 % over the range ± 10 
mm. 

The sensor is isolated from motions of the wind tunnel walls and motions of 
the laboratory by suspending its case in an elastic suspension system (see Fig. 
III-9) with a natural frequency of about 2 Hz. 

F-2 Pitch Sensor 

The pitch angle sensor is a Kearfott type R235 -1-A resolver which is 
excited with a 400 Hz signal. It is linear within 1 % over the range ± 5 
degrees. 

The resolver shaft is connected to the top flexural pivot housing using a 
flexible bellow coupling for protection. The resolver case is held by a bracket 
fixed to the top plunge suspension cross beam. 

F-3 Linear Acceleration Sensors 

Two different types of Piezoelectric Accelerometers, each with a different size 
and sensitivity, made by ENDEVCO are used -- model 2221D shear 
accelerometer weighs 12 gm with sensitivity 17 pC/g, and model 7701-100 
isoshear accelerometer weighs 29 gm with sensitivity 100 pC/g. 
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Both types of accelerometers produce output charges proportional to the 
acceleration inputs. These outputs are converted by a model 2721B charge 
amplifier to voltage readings with a full scale (FS) ± 10 V. 

F-4 Airspeed Sensor 

Airspeed is measured with a United Sensor type PCD-8 KL pitot static 
tube connected to a gage-fluid filled manometer to read dynamic pressure in 
inches of water (a range of 0-8 in. or 0-2000 N/m 2 ). The manometer is cali- 
brated with 0.1 inches graduations (24.9 N/m 2 ). The pitot tube is mounted 
approximately one half-chord ahead of the airfoil leading edge, and a similar dis- 
tance above the test section centerline. This system has an accuracy of 1 %. 
The ambient temperature and pressure are recorded at the wind tunnel for future 
conversion of the dynamic pressure to airspeed. 

G. COMPUTER SYSTEMS 

A DECLAB-23/MNC computer system is used for data acquisition and 
manipulation, monitoring and controlling of external apparatus. Its RT-11 
operating system is a real-time, single-user operating system using FORTRAN 
IV programming language and PDP-11 assembly language. The Real-ll/MNC 
software is designed for use with the RT-11 operating system to support the 
MNC-series hardware modules, which includes the real-time clock, digital 
input/output units, A/D and D/A converters, and a dual- multiplexer. The 
module function diagram is shown in Fig. HI-12. 

The A/D converter is a successive approximation type analog-to-digital 
converter with 12 bits resolution. Its full scale input range is ± 5.12 V bipolar 
with 40 /isec typical data acquisition time. With its internal and the external 
dual-multiplexer it can sample up to 12 channels in quasi-differential mode. 
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There are 4 D/A converter channels with 12 bits resolution and output vol- 
tage ranges ± 2.56 V, ± 5.12 V, ± 10.24 V selected by users. 

A Nicolet 660B dual-channel FFT analyzer is used for system 
identification; it facilitates the measuring of transfer function, damping ratio, 
spectrum distribution, and correlation function, etc. 



Fig. IK- 12 MNC-series Module Function Diagram 


H. SYSTEM PARAMETERS 

The definitions and experimentally measured values of the important physi- 
cal parameters are summarized in Table m-1 along with error estimates. The 
procedures used to determine these values and error sources are detailed in 
Chapter IV and Appendix C. 
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Table HI-1 Definition of Parameters 


Definition of parameters 


Parameter 

Definition 

Experimental 

Value 

Experimental 
Accuracy (%) 
(worst case) 

Units 

/ 

wing span 

0.4191 

0.5 

m 

e 

wing chord 

0.2413 

0.8 

m 

b 

semichord 

■a 

0.8 

m 

a 

distance in semichord 
the EA lies aft of the MC 

-0.3 

0.8 

- 

m 2 

wing mass 

2.1 

9.5* 

kg 


suspension mass 

3.395 

4 

kg 

M t 

total mass 

5.7 

7.2 

kg 

k 2 

mass ratio 

mm 

16.7 

- 

II 

total moment of inertia 
about the EA 

m 

2.4 

Q 


plunge spring rate 

10.38 

1 

kN/m 

K 

pitch spring rate 

43.04 

1 

N-m/rad 


* See explaination in Appendix C. 
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Table HI-1 Definition of Parameters (cont) 


plunge modal frequency 


pitch modal frequency 


plunge damping ratio 


pitch damping ratio 


normalized distance 
the CM lies aft of the EA 


dynamic coupling coefficient 


normalized squared 
radius of gyration 


moment of inertia ratio 


normalized four linkage 
length 


linear motor force constant 


torque motor moment constant 


mass ratio 




0.04814 


.3767 


0.9647 


1.684 


5.295 


0.2232 



3.5 I kg-m 





























































I. APPARATUS PERFORMANCE 

After a thorough investigation of the frequency response of the physical sys- 
tem by using the FFT analyzer, the the overall performance of the system is 
sure to behave very closely to the ideal two DOF typical section shown in Fig. 
III-l. By comparing with the analytical transfer functions obtained in Fig. A-l, 
Section B-l, Appendix A, it can be seen that Figs. Ill- 13 to III-16 truly 
describe them well in terms of the pole-zero locations. However, two complicat- 
ing characteristics exist because of the four-bar linkage used with the torque 
motor and the excitation of differential modes. 

1-1 Unbalanced Force Problem 

Torque is transmitted to .the elastic axis of the airfoil as described in Sec- 
tion E-2. However, the application of a torque generates an unbalanced force in 
the plunge DOF which is 180 degree out of phase with the linear motor force. 
This is explained with the aid of Fig. DI-17. Although this force has been taken 
into account in the system modeling, it may still excite the differential plunge 
mode as shown in Fig. HI- 18. 

1-2 Differential Modes Problem 

The plunge differential mode discussed above which exists due to the 
arrangement of the plunge suspension and asymmetry of the system structure. 
Its existence can be easily removed by adjusting the mass center of the plunge 
DOF; abo, it is not quite observable by the plunge position sensor in our system 
configuration. However, if the amplitude of this differential motion got large 
enough, it could cause binding of the plunge actuator, although this has not 
occurred so far. 
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Magnitude ( dB ) Magnitude (dB) 




Frequency (Hi) 


Fig. m-13 Magnitude of Transfer Function of A /F 


1 



Frequency (Hi) 


Fig. m-14 Magnitude of Transfer Function of a/F 
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Magnitude (dB) Magnitude (dB ) 



Fig. HI-16 Magnitude of Transfer Function of a/ T 
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I 

I 



(It) DIFFERENTIAL 
PITCH. MOTION 



Fig. m- 18 Differential Plunge and Pitch Motions 
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The pitch differential twisting of the wing across its span exists due to the 
- flexibility of the wing spar and the application of torque to only one end of the 
wing as shown in Fig. A-2. It can also be observed in Figs. III-15 and HI-16 at 
about 90 Hz. Since it is much faster (9:1) than the primary system dynamics, 
it can be ignored without too much influence upon our investigation. 

Both differential modes are antisymmetric and thus do not affect the aero- 
dynamics. Since in closed loop studies, it is theoretically possible to drive the 
differential modes unstable, they must be examined closely at that time [26, p. 
83]. 
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Chapter IV 


EXPERIMENTAL METHODS AND RESULTS 


A. INTRODUCTION 

This chapter describes the experimental methods and compares the experi- 
mental results with the theoretical predictions. 

B. MEASUREMENT OF SYSTEM PARAMETERS 

A Nicolet 660B dual-channel FFT (fast Fourier transformation) analyzer 
was used for most of our system parameter identifications. It performs the 
hardware fast Fourier transformations on both input channels, from which a lot 
of functional relationships between these two channels can be obtained. For 
instance, the transfer function between an input and an output of a physical sys- 
tem can be easily obtained by simply plugging its input to channel A and its 
output to channel B. The definition of the functions that the FFT analyzer can 
analyze is given in Table IV- 1. 

With the help of the clamps on both of the plunge and the pitch suspension 
systems, the airfoil can be constrained to move in either single degree of freedom 
(DOF) independently. The FFT analyzer was first used to investigate the 
transfer functions for these single DOF motions. The problems discussed in the 
apparatus performance of Chapter HI were noticed at this stage. The system 
structural parameters were found by looking at those transfer functions with a 
light square aluminum bar installed in place of the airfoil section. 
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Table IV- 1 Definition of Functions for FFT Analyzer 


FUNCTION 

DEFINITION 

instantaneous 
time history 

m 

instantaneous 

S A = F { A(t) } 


T 

FFT 

— lim / A(t)e~ JUJt dt 
T-*oo _j> 

average 

G AA = S A^A 

power spectrum 

= 15/ 

average 
cross spectrum 

Gab = ^b'^a 

transfer 

function 

„ g ab S B -sf A 

n AB = r — 

g aa s A -r A 

impulse 

response 


transmissibility 

r / g bb s b'$b 

m ab — r — '■ 1 

g aa s A -zr A 
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Table IV- 1 Definition of Functions for FFT Analyzer (cont) 


coherent 

G aq' 

cop = n 

G AA 

output power 

— Tab g bb 

coherence 

.2 G AB'(*AB 

Tab /~i s~i 

U AA' U BB 

1 Gj* 

g aa’ g bb 

auto-corr. 


function 

T 

= i im Tt J A (0 A (*+r) dt 

I —►00 Z 1 

cross-corr. 

= r l i Gab } 

function 

1 7 

= lim — [A{t)B{t+T) dt 

i — ♦ 00 Z 1 _ y » 
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Some of the parameters of the airfoil section itself were measured easily out- 
side of the wind tunnel section, such as its mass, dimensions, etc. The available 
data were used to fit the second order clamped single DOF system models in 
order to compute those parameters not easily measured, such as the moment of 
inertia of the airfoil, etc. Redundant data were sometimes available for cross 
checkings. The system parameters measurements and identification are described 
in Appendix C with greater details. 

After all the necessary system parameters were measured or identified, they 
were used in the complete two DOF system model to describe the full system 
behavior. 


C. ACCURACY OF SYSTEM MODELING 

The mathematical system model was derived in Appendix A. It incor- 
porated the unsteady aerodynamics modeling for the 2-dimensional wind tunnel 
airfoil developed by Rock. In his thesis, he found that as H/b decreases, fewer 
poles are needed to be retained to describe accurately the Timman’s modified 
Theodorsen function, i.e. Qj{9,\). If H/b<2, the residue of the first pole 
represents roughly 98% of this function’s impulse response. Consequently, he 
introduced and proved that the approximation 


9+p 


( 4 . 1 ) 


describes Qj{s,\) accurately throughout the entire s-plane for H/b< 2, where 
p is a function of H/b. 
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From Fig. IV- 1, the Nyquist plot of Eq. (4.1) in s/p plane, we can see it 
is a semicircle with center at G+G{zfp-i)f2 and radius G(z/p-l)/2 [Refs. 59- 
60]. The values of G and zjp determines the location and size of this semicircle. 

Playing with the values, we have found a better fit of the theoretical root- 
locus to the experimental one than that found by Rock. Rock chose G = 0.5 
and z/p — 2.0 to preserve the d-c and high frequency gains (1.0 and 0.5 
respectively). Their old and new values are given in Table IV-2. 

A 5th order system model was established to describe the plunge, pitch and 
aerodynamics behavior of the physical system. The resulting modes are named 
by their dominant components for ease of reference. In this model, the unbal- 
anced force due to the applied torque mentioned in Chapter HI-1 was modeled 
into the control distribution matrix G. The damping coefficients were measured 
by the FFT analyzer and added to the modeling. Also the dynamical coupling 
between the plunge and pitch mode was actually computed through the offset of 
the mass center (CM) from the elastic axis (EA) rather than using a rough 
2nd order experimental fit to find the value of the coupling. It is upon the 
center of mass of the wing instead of the center of mass of the whole plunge 
structure that the aerodynamics has influence anyway. These proved to be useful 
to increase the system modeling accuracy. 
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Fig.rV-1 Nyqufat Plot of C?(— ) in - Plane 

*+P P 

Table IV-2 Old and New Values 



old 

new 


value 

value 

p 

0.27 Hz 

0.43 Hz 

G 

0.5 

0.47 

z 

2. 

2 

P 
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A laboratory test instrument, Genrad 2515 structural analyzer, was used to 
check the validity of the system modeling. It gathered a batch of data and then 
used a lattice filter for adaptive processing [61j to fit a specified order linear 
transfer function to the raw data. Figure IV-2 shows a typical output which 
identifies the system eigenvalues directly. 

After repeating the procedure at various wind speeds, those test results could 
be compared with the theoretical predictions by examining the pole-zero loci of 
their transfer functions. Figure IV-3 shows both the experimental and the 
theoretical open-loop system root locus versus various wind speeds. Note that 
they agree within 1 % in frequencies for the worst case before fluttering. This 
result is considered to be better than Rock’s 2 % error. 

It was proven by Rock that a doublet input to the plunge actuator can be 
used to excite both the plunge and the pitch structural modes of the system. 
Actually, a positive voltage (2 volts) was applied to the linear motor for 10 ms 
followed a negative voltage of the same amplitude and duration. This same 
input was used to drive the mathematical model in computer simulations. The 
actual and computer-simulated system responses could be compared in a real 
time manner. Figure IV-4 shows the setup for the real time comparison. 

From Fig. IV-5 to Fig. IV-8, the simulated responses from the computer’s 
digital to analog (D/A) converters were compared with the actual system 
responses. The agreement can be seen excellent for lower wind speeds. However, 
there is a 1 % difference in frequency for the worst case before fluttering shows 
up at a wind speed of about 27 m/sec in Fig. IV-8. 
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Fig. IV-2 Typical Output of the Genrad 2515 Structural Analyzer 
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Fig. IV-5a Actual and Simulated Plunge Responses to 
A Linear Doublet Input at Wind Speed U = 0 m/sec 



Sec 


Fig. IV-6b Overlapped Plot 
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Fig. IV- 7 a Actual and Simulated Accelerometer #3 Responses to 
A Linear Doublet Input at Wind Speed U = 0 m/sec 








D. PARAMETER IDENTIFICATION 


Due to the limitation of our laboratory facilities, only linear and angular 
positions were measured for system monitoring. However, we need all the system 
states to evaluate the estimation error for different sensor locations. Without 
adding other measuring instruments, one easy and possible approach would be to 
use a simple analog differentiating circuit to get their pseudo-rate measurements. 
However, the artificially introduced aerodynamics state was still unknown. 
Finally, we decided to design one more estimator using those position measure- 
ments to construct the pseudo system state for comparison with that obtained 
from the estimator using acceleration measurements. From the discussion in 
Section C about the approximation of using Eq. (4.1) to describe Qj{s,X), it 
seems to be helpful to incorporate some system parameter identification schemes 
to identify those uncertain parameters in order to reduce the system modeling 
error effects. In our case, those key parameters, i.e. p, G, and z/p in Eq. 
(4.1) are good candidates to be identified. Although it was mentioned in Rock’s 
conclusion that the doublet force input did not sufficiently excite the aero- 
dynamic mode for it to be identifiable, an on-line parameter estimation scheme 
through Kalman filtering proposed by Mishne [33] was modified and tried for 
this purpose. The formulations and procedures are followed. 
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First, defining the innovation sequence 


v{k) = {<*) - H%k,a{k- 1)], 

(4.2) 

and the performance index 


J = i u%) W«k), 

(4.3) 

where W is a positive-definite matrix and a(k) b the parameter vector to be 
identified. 

The gradient of the performance index J with respect to the parameter vec- 
tor a(k) is 

m _ w _ Si> T (k) H . 

M j da(k) da(k) * 

(4.4) 

and the approximation to the second order derivative b 


m = 

y ’ da(k) da(k) 

(4.5) 

The on-line application of the Gauss-Newton iteration scheme b 
lowing iteration after each measurement update: 

to make the fol- 

o(jfc) = o(Jfc-l) - At l {k)D(k). 

(4.5) 

After defining the following sensitivity functions 


qiA — dx(k) 

da(k) ’ 

(4.6) 
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(4.8) 

(4.9) 

(4.10) 

then the identification iteration can proceed as follows: 

Given initial values: 

a(0), i(0), 5(0), 

Continuous system modeling: 

i(<) = F\a)x{t) + Gv(t) + Tw(t) 

Discretization and time update of the system state: 

J[*+l,o(*)j = F^a)x\k,a{k)\ + G£a)ti(k ) 

Output predication error: 


and finding the sensitivity matrices 


f (a) = iM 

a[ ’ da(k ) ’ 


FJa) = F d F a (a)T„ 


and 


T, 

Gi, = [/ ^FMri^a, 
0 


f^Ar+l) = s(Ar+l) - Hz[k+l,a{k)] 


Time update of the sensitivity function: 


S(*+l) = F/a)S(t) + + OJ «)«(*) 

Gradient of J: 


D( k+1) = -S T {k+l)H T Wv{k+l) 

Information matrix: 

M(k+l) = S T (k+l)H T WH9(k+l) 


Parameter update: 


o(*+l) = a(k) - M- l (k+l)D{k+l) 

Parameter update of the system state: 

f(A:+l,a(A:+l)] = J{A;+l,a(Ar)] + S(fc+1)[ a(*+l) - a(fc) ] 


Measurement update: 


x(*+l,a(*+l)l = f /- KH)x[k+lXM)] + Ky(k+ 1) 

S(k+1) = ( /- KH)S(k+l) 

If this scheme is used to track the initial condition responses, i.e. v(t) = 0, 
then the computation can be further simplified by not considering the G, G^a), 
and G da (a) matrices. 
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This scheme was checked with computer simulation. To identify the single 
parameter p, it converged within 100 time steps (1 sec) for 20% deviation 
from its nominal value. It converged quickly also for the identification of p and 
G, but had difficulty converging for all p, G, and z/p. 

Unfortunately, the tests performed on the experimental data failed due to 
the existence of the process disturbance and measurement noise since the 
identified parameter values were equivalent to these noise levels. This subject 
will be left for future research. 

E. ESTIMATOR DESIGN 

In this section the estimator design is discussed. The sensor locations are 
evaluated based on the proposed sensor location index. 

E-l Choice of Q w and R v Matrices 

Given an actual design problem, one can often assign a meaningful value to 
R v which is based on the sensor accuracy. The same cannot be said for Q w . 
Physically, Qy, is crudely accounting for both unknown disturbances and imper- 
fections in the plant model. In the design process, we are usually forced to pick 
values of Q w and to settle on an acceptable one after a certain amount of trial 
and error based on the quality of the estimation obtained using gains given by 
specific Q v 's and i?„’s [53]. 

Although we use P as an indicator of relative estimation accuracy, it is 
only an absolute accuracy predicator when we choose the values of Qy, and R v 
based on some knowledge of processes which are approximately white and more 
or less representative of the actual disturbance and noise levels. This makes the 
evaluation of our experiments difficult because the real disturbance and noise are 
not white. Whereas we may be able to find a Q w and R v which gives estimator 
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gains that minimize the estimation error, the error we calculate from this Q w 
and R v may have no relation to measured estimation error. The other difficulty 
we will encounter is that not only the noise levels are different between the two 
different types of accelerometers we used, but also they are different even among 
the same type of accelerometers. Fortunately, the inherent sensor noise level is 
low enough. One easy approach to overcome the above problems in checking 
theory versus experiment is to create a known external disturbance and add a 
known noise to our measurements. 

E-2 Sensor Location Evaluation 

Figure IV-9 shows the setup to evaluate the sensor locations, in which the 
estimator 1 using the position measurements to construct the pseudo full system 
state as the reference and the estimator 2 using the acceleration measurements 
to construct the system state estimate. 

There are only three acceleration measurements available simultaneously due 
to the limited number (three) of the available accelerometer charge amplifiers. 
However, these three measurements can be used altogether to get the estimation 
error of the three sensors case, or separately for single or double sensors cases. 

To simplify the presentation of the test results, both the theoretical and 
experimental values for single and double sensors cases are normalized by their 
corresponding three sensors values. So in the following figures, one will be the 
bottom line which corresponds to the use of all three sensors altogether. In these 
figures, the solid lines are the theoretical predictions which take into account the 
correlation between the process disturbance and measurement noise while the 
dotted lines do not. Generally speaking, the experimental results prove the 
trends predicted by the theory. 
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Fig. IV-9 Setup for the Sensor Location Evaluation 
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F. DISCUSSION OF THE RESULTS 


Figure IV- 10 shows that the experimental data of the first sensor are scat- 
tered about the theoretical value. This may be caused by its improper installa- 
tion. Although this effect is not as significant by looking at the individual esti- 
mation error as in Fig. IV-11. It still can be seen in Fig. IV-12 for the double 
sensor case when one sensor is fixed at the third sensor location. At higher wind 
speeds the error due to this cause is more obvious, see Fig. IV-13. However, the 
result is still good for the other sensors. 


For single disturbance input and single measurement case, the effect of 
changed sensor location amounts to the change of the zero locations of the 
transfer function from the disturbance input to the measurement output. This 
can be seen in Fig. IV-14 which shows the estimator pole locus for increasing 



ratios. The corresponding sensors are numbered from one to six from the 


leading edge to the trailing edge. Also due to the limitation of those trapping 
transfer function zeros, the estimation gain and speed will not be able to increase 
at will. 
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OBSERVABILITY INDEX, trace (WP) 


WIND SPEED 


= 1.00 m/sec 

O-THBORBTICAL VALUES ♦ - EXPERIMENTAL VALUES. NO. 3, 

0* EXPERIMENTAL VALUES. NO. 1. X- EXPERIMENTAL VALUES. NO. ♦, 
A -EXPERIMENTAL VALUES. NO. 3. O - EXPERIMENTAL VALUES. NO. 3. 



- 1.00 - 0.79 - 0.90 - 0.29 0.00 0.29 0.90 0.79 1.00 

SENSOR LOCATION 

Fig. IV- 10 Sensor Location Index vs Different Single Sensor Locations 
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EST. ERROR VS. SINGLE SENSOR LOC. ALONG THE WING CHORD 

WIND SPEED = 1.00 m/sec 

O -THB0RBT1CAL VALUHS X - EXPERIMENTAL VALUES, NO. 3, 

A - EXPERIMENTAL VALUES. NO. 1. © - EXPERIMENTAL VALUES, NO. ♦. 

4* - EXPERIMENTAL VALUES, NO. a. V - EXPERIMENTAL VALUES. NO. S, 



Fig. IV- 1 1 Estimation Error vs Different Single Sensor Locations 
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OBSERVABILITY INDEX, trace (WP) 


WIND SPEED - 1.00 m/aec 

□ -THEORETICAL VALUES + - EXPERIMENTAL VALUES. NO. 3. 

O- EXPERIMENTAL VALUES. NO. 1. X- EXPERIMENTAL VALUES. NO. ♦, 

A -EXPERIMENTAL VALUES. NO. 2. « - EXPERIMENTAL VALUES. NO. S. 

1.3 
1.2 
1.1 
1.0 
0.9 
0.8 
0.7 
0.6 
0.3 
0.4 
0.3 
0.2 
0.1 
0.0 

-1.00 -0.73 -0.30 -0.23 0.00 0.23 0.30 0.73 1.00 

SENSOR LOCATION 

Fig. IV- 12 Sensor Location Index vs Different Doable Sensor Locations 

WITH ONE FIXED AT 0.76 
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OBSERVABILITY INDEX, trace (WP) 


WIND SPEED = 27.19 m/sec 

Q-THIORBTICAL VALUES +-EXPERIMENTAL VALUES, NO. 3, 

O- EXPERIMENTAL VALUES. NO. 1. X- EXPERIMENTAL VALUES. NO.*. 

A -EXPERIMENTAL VALUES, no. z , experimental values, no. s. 

1.2 
1.1 
1.0 
0.9 
0.8 
0.7 
0.8 
0.5 
0.4 
0:3’ 

0 ; 2 ' 

OiE 
0705 

-li00‘ -07735 -0750 -0 23*' 0 00 0.25 0.50 0.75 100 

SENSOR' LOCATION 

Fig. IV-137 SensorLocationlndexvs Different Single Sensor Locations 
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Chapter V 


SUMMARY AND RECOMMENDATIONS 


A. SUMMARY 

The problem of placement and number of sensors has been investigated. A 
sensor location index based on the weighted sum of the mean square estimation 
errors is proposed and has been experimentally verified. Although there have 
been many proposals regarding the choice of the criterion, the sensor location 
index features easy computation, having direct physical interpretation, and show- 
ing relative improvement for an increased number of sensors. Also this new cri- 
terion makes the gradient formulation easy to handle and compute. A computer 
program for the systematic optimal sensor location search is developed based on 
the gradient search optimization technique. 

Since the use of accelerometers as measuring devices is so common in most 
control applications, the control design computer program OPTSYS has been 
modified to accept a more general input form including a rate measurement and 
which also takes into account the correlation between the process disturbance 
and measurement noise introduced by those accelerometers. Similarly, an itera- 
tion scheme doing this for the discrete case is also developed to facilitate the esti- 
mator implementation on digital computers. 

An experimental study for choosing the sensor location has been conducted 
on an aeroelastic system. It consists of a NACA 0015 typical section wing with 
six accelerometers installed inside along the wing chord as the esti- 
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mator measuring instruments, an existing wind tunnel section, and some other 
accompanying experimental devices. Fabrication techniques are developed and 
apparatus built for experimental studies. 

A sensor location for which the system is unobservable has been verified to 
be the center of percussion of the rigid two degree of freedom typical section 
wing. System modeling which includes the unsteady aerodynamics model is 
improved and experimentally verified. An on-line parameter identification 
through Kalman filtering has been formulated and applied to aerodynamic 
parameter identification. It is successful in simulations when no system distur- 
bance and measurement noise exists, but fails for the actual system test due to 
the system noises. 

Sensor location choice is evaluated experimentally based on the proposed 
sensor location index. It proves the trend predicted by the theory. This experi- 
mental evaluation on an aeroelastic system is significant since it is believed to be 
the first for feasibility studies. 

After investigating the double sensor example which results in collocated 
sensors at the trailing edge we can see the negative side of the sensor location 
index. This result shows that it may be too dependent on the theoretical model- 
ing of the system and it suggests that a more realistic index should also take into 
account the parameter uncertainties of the system for sensor location selection. 
Specifically, one would like to separate the sensors if there was no model 
knowledge and depend only on the kinematic information available from the two 
sensors for the two mechanical degrees of freedom. The results also show the 
need for accurate models for the process disturbance and measurement noise 
when making critical decisions. 


- 109 - 



B. RECOMMENDATIONS FOR FUTURE RESEARCH 

1) Newton-Raphson technique can be applied to the sensor location search pro- 
gram for faster searching speed. Better schemes should be tried to detect 
the local minimums. 

2) Extend the same idea to the optimal actuator location search through the 
duality between the estimator and controller design. 

3) Estimation and sensor location analyses in the presence of system parameter 
uncertainties and modeling errors should be investigated. This should also 
include the exploration of a more realistic criterion which takes greater 
account of the kinematic information of the system as model uncertainties 
increase by considering a frequency-band weighted SLI. 

4) Further investigation of system parameter identification with the current 
setup is suggested using optimal input design. 

5) Adaptive estimation with combined on-line system parameter identification 
and state estimation should be developed. 

6) Apply these results to the airfoil section with a trailing-edge flap (built by 
Stoltz) for improved estimation and feedback control of flutter. 
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Appendix A 


t 

MATHEMATICAL MODELING 


In this Appendix, Section A is given for completeness of the derivations of 
equations of motion, although it is well known in many textbooks [46]. Section 
B includes some refinements to differentiate the mass of the airfoil and the mass 
of the whole plunge system. Section C and Tables A-l - A-5 are summarized 
from Rock’s thesis [25] for ease of reference and they are corrected for typo- 
graphical errors. 

A. EXTENDED HAMILTON PRINCIPLE 

In general, for a system of N particles in equilibrium, the sum of the vir- 
tual works over all particles must be zero, or 

6W = £ SW t = E Hi**.* = 0- (A-l) 

n=l n=l 

where R,- is the resultant force vector acting on a particle at position r,- [46]. 

After eliminating all the constraint forces and the internal forces among 
rigid bodies, we obtain 

E F .- • 6r i = °> (A-2) 

n=l 

where F,- is the resultant vector of externally applied forces at r,-. 
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Equation (A-2) is the expression of the virtual work principle. It can be 
stated as: If a system of forces is In equilibrium, the work done by all the 
externally applied forces through virtual displacements compatible with 
the constraints of the system is zero. 

By incorporating D'Alembert's principle, which states that the resultant 
force is in equilibrium with the inertia force, one can extend the principle 
of virtual work to cover the dynamic case following the same reasoning for Eq. 
(A-2), i.e., 

E (F,- - p,) <5r- = XI (F,-- m t r { yST { = 0. (A-3) 

n= 1 n= 1 


Note that 


E F,- • 6r; = SW t 

n«=»l 


(A-4) 


if, * 6r i = fai ’ hi) (r,- • r,), 


(A-5) 


and 


XI w,r- • 6r { = XI m i 4l(*i • (5r .) - s E ' *•*) 


n=l 


n=l 

= E m ,- • &\) - *7- 


(A-6) 


Introducing Eq. (A-4) and Eq. (A-6) into Eq. (A-3), one obtains 


6T+6W 


N d 

E m i 4-A T i • <5r,). 

n=l 


(A-7) 


Integrating Eq. (A-7) between <j and f 2 , the result is 
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(A-8) 

(A-9) 


Equation (A-8) is referred as extended Hamilton principle by letting W 
consist of both conservative and nonconservative works. Now Eq. (A-8) 
becomes 



The expressions T and V can be thought of in the following forms (in our 

special case) 


T = T(t), 

(A-12) 

V = V (r), 

(A-13) 

where r is a generalized coordinate vector; and noting that 


6T = 6t, 

dr 

(A-I4a) 

XT/ W , 

SV = — St, 

dr 

(A- 14b) 
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SWnc = Q« • 


(A- 14c) 


where Q ne is a generalized nonconservative force vector. 

After substituting Eqs. (A-12)-(A-14) into Eq. (A-10), integrating by parts 
and using the assumption of Eq. (A-9), also noting that the resulting equation 
equals to zero must hold for arbitrary 6t, we obtain the general form of the 

Lagrange’s equation 

i'f > + f = < A -> 5 > 

B. EQUATIONS OP MOTION 

B-l Two DOF (degree of freedom) Model 

Figure A-l shows the physical system and its modeling representation. In 
this case, 


,T _ 

[ * ° ]■ 

(A-16) 

T = 

jMjif + A 2 + S a ah, 

(A-17) 

V = 

+ \k„ 0 \ 

(A-18) 


The virtual work term arises from the external energy input into the struc- 
tural system by the torque T a about the elastic axis (EA), the force F to the 
EA, and the unsteady pressure difference between the upper and lower wing sur- 
faces [18, p. 39j. The unsteady pressure distribution can be defined as 

Ap (x,t) = p (z,y=0 *,t) - p (x,y=0V), (A-19) 
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Fig. A-l Two Degree of Freedom Modeling: Plunge and Pitch 
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where 


p (x,y=0 + ,t) = static pressure at airfoil upper surface, 
p ( x,y=0~,t ) — static pressure at airfoil lower surface. 

So the first variation of the nonconservative virtual work can be expressed as 


where 

and 


Substituting Eqs. (A-16)-(A-18) and (A-20) into Eq. (A-15), and adding 
the structural damping terms, the complete aeroelastic equations of motion 
become 

(A-23) 
(A-24) 

(A-25) 
( A-26) 

(A-27) 


T 

M-jh + S a d + C h h + K h h = F--y-L (f), 


S a h + I Q 'a + C a a + K a a = T a + M a (<), 


where 


S a — M2X al 

b 

L (f) = - / Ap (z,t)dx, 

-b 

b 

M a (t) = / Ap {x,t){x-a)dx. 
-b 


6W ne = - / Ap (x,t)6ydx + (F —)Sh + T Q 6a, (A-20) 

-b 1 

y = - { h + (x-a) a ], (A-21) 

6y = = -[6h + (x-a)6a ]. (A-22) 
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B-2 Three DOF Mode! 


Figure A-2 shows the same physical system but one more DOF is added to 
it’s modeling representation. In this case, 

r T — ^ h a 2 a\ j, (A- 28) 

T = jM r h 2 + + S a a 2 k, (A-29) 

^ + ~<4 + fa* ~ «,) 2 + ~ o ,\. 

(A-30) 


And similarly, 


6W. 


ne 


0 rp 

-l Ap (x,t)6ydx + (F- -f-)6h + TJ a lf 
-b 1 


(A-31) 


where 


y = - [ h + (z-a)a 2 ]. (A-32) 


and 


6y 





0a 2 ~ - [ 6h + (x-a)5a 2 ]. 


(A-33) 


After substituting Eqs. (A-28)-(A-31) into Eq. (A-15) and adding the 
structural damping terms, we obtain the complete aeroelastic equations of motion 

T 

Mjh + S a d 2 + C k h + K k h = F j~ ~ L (t), (A-34) 
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(A-35) 


■S’* A + 4 a a 2 + ^0^2 + “2“°2 + Ke( a 2 ~ Q l) ~ M* (0> 


K> 


/a^l + + "T" a l + ^f( a l “ °2) — T'a- 


(A-36) 


This model is presented for reference only and not used for analysis. 
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Fig. A-2 Three Degree of Freedom Modeling: Plunge, Pitch and Twist 
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C. UNSTEADY AERODYNAMICS MODELING 
If an experimental airfoil with finite thickness is placed in a flow field 

bounded by tunnel walls, we must consider the blockage and lift interference 

\ 

effects. 

\ 

C-l Blockage Corrections 

Blockage corrections account for the reduction in available cross-sectional 
area when a body of finite thickness is placed in a bounded flow field [62, p. 
300]. These corrections may be divided into solid blockage term and wake block- 
age term for a finite thickness airfoil and its wake, respectively, see Fig. A-3. 
They amount to a correction on the airspeed, i.e. 

U = ( 1 + t lb + e ub ) U ve , (A-37) 

where U ue is the uncorrected free stream flow velocity. The results for our 
NACA 0015 airfoil with a b/H ratio of 0.5278 is a 1.9% increase in the 
airspeed. 


//////////// ///////////////////////// /////////////////// 
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Fig. A-3 Blockage Effects 
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C-2 Lift Interference Effect 

This effect is caused by the wind tunnel walls surrounding the airfoil. 
Assume an infinitely thin airfoil performing oscillations with an infinitely small 
amplitude. Its magnitude is prescribed in every point of the wing chord instead 
of the airfoil contour for simplification. The wind tunnel is considered to be 
infinitely long with constant cross section so that it does not influence the undis- 
turbed flow. The airfoil spans the height of the tunnel section equidistant from 
and parallel to the side walls. The flow field can be considered as two dimen- 
sional and incompressible by assuming invicid, low speed air flow. Under these 
conditions, Timman solves the problem with an exact solution [20] using the 
velocity potential functions and conformal transformations in a manner analogous 
to Theodorsen’s treatment of an airfoil in free (unconstrained) flight [19]. 

Using the coordinate systems in Fig. A-4 and the two transformations 


where 


X+ i Y = 


sinh nz/2H 


sinh xb/2H 


Z = cn (?,**)> 


z = x + i y, 

? = f + * n, 

k * = tanh nb/2H, 


(A-38) 

(A-39) 


and cn is a Jacobian elliptic cosine function, Timman maps the region in the 
2 -plane bounded by the airfoil and the wind tunnel walls first into a cut 2^-plane 
and then into a rectangle in the f-plane. These two successive transformations 
reduce the problem to the one which he can solve for the disturbance velocity 
potential 4> = + $ 2 . 
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<t>i is the velocity potential for the noncirculatory flow field and satisfies the 
boundary conditions in the 2 -plane 


and 


a<*>! 

dn 


i^x)e Mt : y = 0, - b < x < b; 


dn 


0 : y=±H, 


(A-40) 


where n is in the normal direction to the wing chord or tunnel walls. The term 
v(x) is a known function which describes the motion and shape of the airfoil. <J> 2 
is the velocity potential for the circulatory flow field resulting from the vortices 
along the wake, and subjects to the boundary conditions in the 2 -plane 

dn 

After applying the Helmholtz condition of persistence of vorticity and the 
Kutta condition of finite velocity at the trailing edge, $ can be solved 
analytically. The pressure distribution over the airfoil can thus be found from 
the unsteady Bernoulli’s equation in the linearized form [63, p. 226] 


0 : y = 0, - b < x < b] and y — ± H. (A-41) 


Ap (x,f) = " P ( Jj ( *+ ~ *- ) + u ( «+ - «- ) ]. ( A ’ 42 ) 


where + and - represent the upper and lower surface of the airfoil, respectively. 
Finally, the lift and moment about EA can be obtained by integrating the pres- 
sure distribution over the airfoil using Eqs. (A-26)-(A-27). 

Rock proved that the expressions of the Fourier transform of the lift and 
moment about EA of Timman’s results could be written as [25, p. 27] 
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(A- 43 a) 


L (u>) = a^irpb 2 [ (kj) 2 A + U(ij)a - ba(iuj) 2 at ] 

+ a 2 t l'KpUbQ T (k ,\ ) { ( iu)h + Ua + 6(o3-a)(iw)a ], 

and 

M a (w) = fljjrpfc 2 [ ba{uj) 2 h - 6 2 (a 4 +a 2 )(to;) 2 a - Ub(a$-a)( iui)a ] 

+ 022irpUb 2 (a+a 5 )Q T (k ,\ ) ( ( iu)h + Ua + 6(a3-a)(«u;)a J, (A- 43 b) 

where Q T {k,\) is the modified Theodorsen function. 

The technique of simply replacing ik with Laplace variable a is used to 
extend the modified Theodorsen function to include the general motions of 
an airfoil. 

Definitions of Q T and coefficients a,- are repeated here in Table A-l for 
easy reference. It is also accompanied with free flight values, which are the limit- 
ing values as wall separation approaches infinity, i.e. H/b -*■ 00. Values of 
coefficients a { for different H/b ratios are computed by Rock as in Tables A-2 
- A- 3 . 

Rock also shows the approximation 

~ G ]^ = g + W’ ,a - 441 

describes < 3 r(a,X) accurately throughout the entire a-plane for H/b < 2. 
After incorporating this approximation into Eq. (A- 43 ) and knowing that 

03 (s,X)^ r («,X) = 07 - ag[ 1 - T\3,\) ], (A-45) 

combining coefficients and using new defined coefficients in Table A- 4 , the 
expressions of lift and moment about EA may now be written as 
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L (a) = a x irpb* [ £h + Usa - a^a ] 


4- afixpUb* | [ a n sh 4- a n Ua 4- a 12 sa ] 


+ C(a,X) [ a 10 sh 4- a lQ Ua 4- a 13 aa ] 


(A-46a) 


and 


M a (a) = a x irpb A [ a£h - U[a b -d) 8 Ct - (a 4 +a 2 )a 2 a ] 

+ 022irpUb A (a+a b ) | [ a n sh + a n Ua 4- a 12 aor ] 

4- Cfa.X) [ a l 0 8 h+ a 10 Ua + a 13 aa ] j, 


where 


5(»,x) = HiM+J. = d ipdl 

2 (*+?) 


(A-46b) 


(A-47) 


Introducing an augmented aerodynamic state y as 


V = 


(a+p) 


( a l0 sh 4- a l(i Ua 4- a 13 aa ], 


(A-48) 


and combining Eqs. (A-23)-(A-24) and (A-46)-(A-48) yields 


d x s 2 h + d%sh + d$h 4- d^a + d$sa 4- d 10 a = 


F To 


A/j'ft M T bl 


- d n% 


Aa 2 ® + fs sa + k a + U^ h + Ao*^ = -f- + AiF. 


(A-49a) 

(A-49b) 


where coefficients d { and A are defined in Table A-4. After solving Eq. (A-49) 
for s 2 ^ and a 2 ®, and replacing the Laplace operator a by the time derivative, 
we get Eq. (A-50) with coefficients </,• defined in Table A-5. 
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Table A-l Definition of Coefficients 


























Table A-l Definition of Coefficients (cont) 


























Table A-l Definition of Coefficients (cont) 


















Table A-l Definition of Coefficients (cont) 
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Table A-2 a; at Various H/b 


a,- at Various H/b 

H/b 

OO 

10 

2 

1.745 

1 

0.5 



1.0021 

1.0488 

1.0631 

1.1725 

1.5191 

<*2 


1.0041 

1.0952 

1.1224 

1.3209 

1.8971 

«s(0) 


0.4990 

0.4788 

# 

0.4438 

0.4004 

«4 

m 

m 

0.1168 

0.1148 

0.1023 

0.0827 



m 


0.4670 

m 

0.3751 

II 


0.4990 

0.4788 

0.4736 

m 

0.4004 

II 

<3° 

0.2500 

0.2478 

0.2067 


0.1429 

0.0767 

<3=L 

II 



1.0000 





* See Table A-3. 


- 130 - 
























































Table A-3 Dependence of Timman’s a 3 (k,X) on k and H/b 


Dependence of Timman’s 03 (Ar,X) * on hand H/b 


k H/b 


0.4990 
+ 0i 

0.4788 
+ 0i 

0.4438 
+ 0i 

0.4004 
+ 0i 

0.4994 
+ 0.0008i 

0.4810 
+ 0.01 16i 

0.4470 
+ 0.0229i 

0.4040 
+ 0.0297i 

0.4999 
+ 0.0012i 

0.4867 
+ 0.0210i 

0.4559 
+ 0.0420i 

0.4143 
+ 0.0569i 

0.5003 
+ 0.0013i 

0.4941 
+ 0.0273i 

0.4687 
+ 0.0575i 

0.4296 
+ 0.0798i 

0.5006 
+ 0.0013i 

0.5017 
+ 0.0307i 

0.4832 
+ 0.0683i 

0.4479 
+ 0.0974i 

0.5009 
+ 0.0013i 

0.5085 
+ 0.0319i 

0.4978 
+ 0.0749i 

0.4674 
+ 0.1098i 

0.5011 
+ 0.0013i 

0.5143 
+ 0.0318i 

0.5114 
+ 0.0781i 

0.4866 
+ 0.1177i 

0.5014 
+ 0.0012i 

0.5228 
+ 0.0298i 

0.5339 
+ 0.0782i 

0.5209 
+ 0.1235i 

0.5016 
+ O.OOlli 

0.5283 
+ 0.0271i 

0.5503 
+ 0.0739i 

0.5482 
+ 0. 12111 


* This parameter value is not used, it is for reference only. 
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Table A-4 Definition of Coefficients 


Definition of Coefficients, </,• and k 

dx 

1 + aju 

<h 

2^ /n 


2a 2 a u U/n + 2 $ h u h 

d% 

<k + djd a 10 

<k 

4 

d» 

d b + djd a 13 

d* 

K 2 x q - a x a/n 

dio 

<k + <hd a l0 U 

d, s 

a x U/n + 2 o 2 a l2 U/n 

dn 

(z - P ) 

^6 

20^0^1? I\i 



a 

a 1 (a 4 + a 2 ) 



k 

a i( a 5 - a)U 2 <t 2 (a + a 5 )a l2 U 
K 2 firl K 2 nri ** " 



k 

2(u z {a + a b )a n T? 
or; 

KiV’i 



k 

xJd - a i a /K 2 ni~ t 

k 

h~h^ fl i3 

k 

-2a 2 {a + a b )a u U/K 2 nil 

h 

k~ h(* fl io U 

k 

not used 

Ao 

k~ fad fl io 

h 

2(^0 + a b )U/K 2 nri 

hi 

h<* (* ~ P ) 
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Table A-5 Definition of Coefficients 


Definition of Coefficients, g ,• 

9i = 4 bM “ / 10//4 

II 

to> 

II 

92 = 9s/ 97 

93 — dj d x - / 1//4 

93 = 9n/ 97 

<74 — ^ 9/^1 _ 

94 =: 9io/ 97 

% = <W <*i ~ h/h 

95 = "912/ 97 

% = ~d xx /d x - / 11 /A 

_<§? 

II 

<§? 

< 

I 

■w" 

II 

^S? 

C? 

II 

£ 

98 — ds/d 4 - / 10//1 

98 — 95 / 93 

9fl = d^l d 4 

c? 

II 

<§ 

< 

-5? 

I 

II 

O 

•■n 

©• 

_£> 

1 

II 

O 

*5 

9n — <W d A - / 8 //| 

9n — ~ fl io 

9l2 = ~d\i! d\ - A 1//1 

912 =- fl io^ 

913 — _fl 13 

914 = P 

915 = V <*497 

916 = l//l97 

917 = 1/ <*193 

918 = I //493 
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Table A- 8 Derivative of Parameters 


Derivative of parameters 
w.r.t. G, p, and 

d { and /,• 
z 

d<k 

A 

dd n 

<h( z- p) 

dG 

<h 

dG 

a 10 


^7°13 

dd n 

djG 

dG 

fl 10 

dp 

a io 

dd l0 


dd u 

djG 

dG 

drjU 

dz 

a 10 

df s 

h a i3 

dfn 

hi*- P ) 

dG 

fl 10 

dG 

«10 

dd$ 

-f 7 U 

Vu 

f 7 G 

dG 

dp 

a io 

dfio 

r 

dfn 

f 7 G 

dG 

h 

dz 

a io 
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Table A- 6 Derivative of Parameters (cont) 



Derivative of parameters <7,- and 
w.r.t. G, p, and z 


f 7 

■ • # 


+ h\ JJis 


d\ / 4 J flio 

4+41 v 


^ Al z - p 

d\ fa I a 10 


^7 /7' 

■ * • 


. /7I a i3 


<h | / 7 | G 
d\ / 4 J flio 




^7/7 
« 1 # 








dj fj 

I ' # 


^ 7/7 
> * * 


— a. !l 

dx U\ 

d, fx\ 

— A. — 

a * a 


— a. !l 

a * / 


^7 A 
• 1 • 


a 13 1 
°10 % 


2 - p 1 

fl 10 93 


°13 1 

a io 97 


z - p 1 

a 10 97 



4+4l-£- 


/J fl io 
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Appendix B 


ACTUATOR AND SENSOR CALIBRATIONS 


Characteristics of sensors and actuators enter into the determination of sys- 
tem parameters (see Appendix C) and are important in predicting system 
dynamic behavior. Accurate calibration is therefore necessary. 

Since the same test section built by Rock was used for our experiments with 
a new airfoil, most of the calibration procedures were similar to those used by 
him, however, with some different approaches. While every effort was made to 
reproduce or cross check with Rock’s results, the calibrations were carried out 
independently with different tools and equipment. 

A. ACTUATORS 


A-l Static Response 

The static calibration procedures of the plunge and pitch actuators are out- 
lined in Fig. B-l. A PDP-11 computer system was used to implement the digi- 
tal controller with the actuator as the controlled element. Known disturbance 
loads (weights hung in the gravity field) were balanced by the action of this 
controller. A 60 Hz dither signal was added to eliminate the effects of stiction. 
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Fig. B-l Experimental Configuration for 
Static Calibration of Actuators 
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A-2 Dynamic Response 

A Nicolet 660B dual-channel FFT analyzer was used with a Tektronix 
P6042 current probe to find the frequency responses of the actuators and the 
current-drive amplifier. 

A Hewlett-Packard 3722A noise generator was used to generate an infinite 
sequence, 150 Hz bandwidth (system dynamics was about 10 Hz) noise signal 
to drive the amplifier. First the transfer functions between the current outputs 
from the amplifier and the voltage inputs to the amplifier were obtained using a 
pure resistive load. It verified that the amplifier worked as a constant gain 
amplifier (K Aplunge = 0.590 A/V, K Apitch = 0.270 A/V) with zero phase 
shift within 0-100 Hz. Similar transfer functions with actuators connected to the 
amplifier were obtained, which showed a small change in amplifier gains 
(K A plunge = 0.630 A/V, = 0.273 A/V) because of the back electro- 

magnetic force (EMF) from the actuators. 

Then the transfer functions between the position sensor voltage outputs and 
the actual actuator current inputs were obtained. These transfer functions 
showed the responses of the system also (see Fig. B-2 and B-3). Above 10 Hz, 
the actuators should behave as forces (torques) acting on pure masses (moment 
of inertias) with minus 40 dB/decade magnitude slope and 180 degree phase 
lag. However, there were significant phase shift and change in magnitude slope. 
Specifically, a big magnitude drop and phase shift at 63.5 Hz in the plunge 
motor’s response and 45 degree phase shift at 60 Hz in the pitch motor’s 
response. They were caused partly by the associated system dynamics and those 
high frequency resonances of the actuators, as well as their bandwidth limita- 
tions. 
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Phase (degree) Magnitude (dB) 



Fig. B-2 Frequency Response of h / F 
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I'hase (dB) 



Frequency (Hi) 



Fig. B-3 Frequency Response of afT 
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B. SENSORS 


B-l Static Response 

The plunge displacement sensor, LVDT, was calibrated by clamping its 
core in the chuck of a milling machine (Societe Genevoise, SIP, Switzerland) 
and fixing its case to the movable table. Table motion could be read to 
2.54 *10" 3 mm (10 -4 in.). Figure B-4 shows the voltage output of the LVDT 
versus the known table displacement using linear least square fit to the raw data. 
Its gain constant is 15.031 V/in. (627.2 V/m) which is accurate to 1 % over 
the ±10 mm limit on plunge displacement. 

The pitch displacement sensor, resolver, was calibrated with a dividing head 
on a mechanical lathe (Leitz, Wetzlar, West Germany). It could be rotated 
accurately down to 1 arc second ( 2.78 ♦lO -4 degree). The case of the resolver 
was fixed while its shaft rotated with the dividing head. Figure B-5 shows the 
resolver voltage output versus known dividing head rotation using linear least 
square fit to the raw data. Its gain constant is 52.13 V/rad which is accurate to 
1 % over the range ±5 degrees. 

B-2 Dynamic Response 

No tests for either sensor due to the difficulty to perform with available 
equipments. However, we anticipate no effects would occur below the 400 Hz 
excitation frequency. 
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Fig. B-4 LVDT Calibration Curve 
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Fig. B-5 Resolver Calibration Curve 
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Appendix C 


MEASUREMENT OF SYSTEM PARAMETERS 


Brief discussions of the procedures used in determining system parameter 
values are presented with estimates of accuracy (worst case). 

Figure C-l shows the dimensional parameter definitions with exaggerating 
variations for clarity. 


SPAR 



Fig. C-l Plane View of Wing Exaggerating Dimensional Imperfections 


- 144 - 



t = 0.4191 m: wing span. 


Directly measured with an accuracy of 2 mm (0.5%). 

c = 0.2413 m (6 = 0.12065 m): wing chord (semichord). 

Template was numerically machined with an accuracy of 0.0254 mm 
(10 -3 in.) (0.01%). Form section cut with hot wire has imperfections in 
the leading and trailing edge may result in chord variations of 2 mm 
(0.8%) along the span. 

a = - 0.3: distance in semichords the elastic axis (EA) lies after the mid- 
chord (MC). 

Variations in a are of the same order as variations in c (0.8%). 

M m = 0.3546 kg, 0.8461 kg/unit span (m): Al spar mass. 

Directly measured using a digital scale with an accuracy of 10 -5 kg 
(0.003%). 

M 2 = 2.1 kg, 5.011 kg/unit span (m): total mass of the wing section. 

Directly measured as M Al . It included the mass of plexiglass endplates 
(0.5473 kg, 26.1% of Af 2 ), fixing screws, accelerometer connecting 
wires, etc. The way the connecting wires being handled may result in 
0.2 kg (9.5%) variations in Jl/ 2 . 

u}^ Al = 8.375 Hz, 52.62 rad/sec; ^ Ai = 0.0152: uncoupled plunge natural 
frequency and damping ratio with wing section replaced by Al spar only. 
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A Nicolet 660B dual-channel FFT analyzer was used to find this 
plunge mode frequency response with pitch mode clamped. Resonant 
peak can be determined accurately within its resolution range (0.125 
Hz, 1.5%). Damping ratio was computed by the analyzer using the 
relationship f = 1/2 Q, where Q was the amplification factor 

2 = /o/(/ 2 -/i)- 

u a,Ai — 21.25 Hz, 133.5 rad/sec; $ a A i = 0.0148: uncoupled pitch natural 
frequency and damping ratio with wing section replaced by Al spar only. 

Same method as above, which is accurate to 0.125 Hz (0.6%), except 
the plunge mode was clamped. 

u h,xcin 9 = 6 75 Hz, 42.41 rad/sec; ? ^wing = 0.00111: uncoupled plunge 
natural frequency and damping ratio with wing section installed. 

Same method as above, which is accurate to 0.05 Hz (0.7%). 

Voting — 9.70 Hz, 60.95 rad/sec; S a ,wing — 0.00981: uncoupled pitch 
natural frequency and damping ratio with wing section installed. 

Same method as above, which is accurate to 0.05 Hz (0.7%). 

u k = 6.65 Hz, 41.78 rad/sec; = 0.00448: coupled plunge modal fre- 
quency and damping ratio with wing section installed. 

Same method as above, which is accurate to 0.05 Hz (0.7%). 

u a = 10.05 Hz, 63.15 rad/sec; = 0.00932: coupled pitch modal fre- 
quency and damping ratio with wing section installed. 

Same method as above, which is accurate to 0.05 Hz (0.7%). 


- 146 - 



K h = 10.38 kN/m, 24.78 (kN/m)/unit span (m): plunge spring rate. 

Measuring the displacement of the plunge motion due to a known 
applied load (weight) yielded a direct measure of the plunge spring 
characteristics (Fig. 02). Accuracy depends on the calibration of the 
plunge displacement sensor (1%) and the weight (0.01%), totally 
equals 1%. 

K a — 43.04 N-m/rad, 102.7 (N-m/rad)/unit span (m): pitch spring rate. 

Measuring the angular displacement of the pitch motion due to a known 
applied torque (weight hung on one end of the four-bar linkage) yielded 
a direct measure of the pitch spring characteristics (Fig. 03). The 
angle measured by the resolver was actually the twisted wing angle, 
however, this fact was ignored in this measurement. Accuracy depends 
on the pitch displacement sensor (1%). 

M TAl = 3.750 kg, 8.947 kg/unit span (m): total apparent mass with Al 
spar installed only. 

The apparent mass was found from and as 

Mt,ai - -r-. (C-l) 

< 4 , Al 

It included the mass of the Al spar and the suspension systems. Its 
accuracy depends on K * (1%) and (1.5%), and is about 1 + 

1.5 * 2 = 4%. 
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Fig* C-2 Linear Spring Constant Calibration Curve 
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F!g. C-3 Angular Spring Constant Calibration Curve 
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M x — 3.305 kg, 8.100 kg/unit span (m): total apparent mass of the suspen- 
sion systems. 


A/j = M T Al - M Ab 


( 02 ) 


which is accurate to 4%. 


/ 0l = 1.502*10 -3 kg-m 2 , 3.584*10 -3 kg-m 2 /unit span (m): apparent 

moment of inertia of the pitch suspension system about the EA. 

This apparent moment of inertia about the EA was found from K a 
and as 


i 



( 03 ) 


Since the calculated moment of inertia of the Al spar was only 1% of 
I ai , it could be ignored. I Qi is accurate to 1 + 0.6 * 2 =s 2.2%. 

I a — I at — 1.151*10 -2 kg-m 2 , 2.746*10" 2 kg-m 2 /unit span (m): total 
apparent moment of inertia with wing section installed about the EA. 

4 r = -r 2 -’ (C* 4 ) 

^ Q,mng 

which is accurate to 1 + 0.7 * 2 = 2.4%. 

4 a = 1.001 *10“ 2 kg-m 2 , 2.388*10 -2 kg-m 2 /unit span (m): total apparent 
moment of inertia of the wing section and endplates about the EA. 
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4j 4r 4i» 


(C-5) 


which is accurate to 2.4%. 

M t = 5.7 kg, 13.60 kg/unit span (m): total apparent mass of the wing sec- 
tion and the suspension systems. 

Mf — M\ + Mty (C*6) 


which is accurate to 7.2%. 

7^ = 0.3767: squared radius of gyration in semichords squared about the 
EA. 



4 T 


(07) 


which is accurate to 2.4 + 7.2 + 0.8 * 2 = 11.2%. 

x a — 0.19: distance in semichords the. center of mass (CM) lies after the 
EA. 

The measurement of x a is based on its physical definition shown in 
Fig. C-4. Oscillation frequency was measured with a stop watch accu- 
rate to 0.1 sec (0.3%) for 30 cycles (27.5 sec). 



(08) 


which is accurate to 2.4 + 0.3 * 2 4- 7.2 4- 0.8 = 11%. 
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I Q o t + M 2 gx a a = 0, 
w 2 / a = M 2 gz a , 

M 2 g- 

Fig. C-4 Calibration of x a 

T= 1.684: length of the four-bar linkage in semichords. 

Directly measured with an accuracy of 2 mm ( 1 + 0.8 

Ko — 0.3684: mass ratio. 



K = 0.9847: moment of inertia ratio. 



= 1 . 8 %). 


(C-9) 


(C-10) 
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K e = 421.1 N-m/rad: spar coupling spring rate. 


Computed value for 3 DOF (degree of freedom) modeling. 
Kp = 5.295 N/V: linear motor force constant. 

K t — 0.2232 N-m/V: torque motor moment constant. 


o; ai = 91.50 Hz, 574.9 rad/sec; = 0.0443: differential pitch modal fre- 
quency and damping ratio. 

Directly measured by FFT analyzer which is accurate to 0.25 Hz 
(0.3%). 

<jj eu = 6710 Hz, 42,166 rad/sec: coupling frequency from wing to four-bar 
linkage. 



(C-ll) 


w Ci2 = 44,629 Hz, 280,414 rad/sec: coupling frequency from four-bar link- 
age to wing. 



H = 246.0: mass ratio. 


M T /r 

wpl? 


(C-13) 
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The density of air in the tunnel p is nominally 1.219 kg/m 3 , which is 
known within 1%. p is thus accurate to 7.2 + 0.5 + 1 + 0.8*2 = 
10.3%. 

S Q = 0.04814 kg-m, 0.1149 kg-m/unit span (m): dynamic coupling 

coefficient. 

S a = M 2 x Q = M 2 x a b y (C-14) 

which is accurate to 3 + 0.5 = 3.5%. 
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Appendix D 


TABLE OF GRADIENTS 

1. 4 tr( - ] * L 

2. trf A X ) « A # 

•J. ^t r (AX']«A 

4. ^ tr(A X B] • A'B' 

5. *^^fAX»B)»BA 

6 . ^trlAX] « A 

7. ^#tr{AX'] - A' 

H. "§x ,tT l - * £ ] * B A 

9. *^»tr(A X'B] « A'B' 

10 . " 5 ? «r[X X] >2X* 

11. ^ttfXX'l »2X 
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12. 


ABLE OF GRADIENTS (coat) 


| f «rIX n j«n<X n *‘)' 


~ ir j A X" )« [ ^ X 1 A X 0 ' 1 ' 1 ' 


' i =0 

14. Jg rr|A XfiX] * A'X'B* + B»X'A* 

15. -jfa rrf A X B X' ] * A'X B' t AXB 

D , X t X 

16. tr l c J * c 

17. tr f X* 1 ] * -(X‘V V * * (X* 2 )' 

18. tr(A X“ *B 1 * - ( X " 1 B A X " 1 ) * 

19. -^det[X] ■ (der [X ]) (X* V 

20. -2- log ctet f X J *<X‘ V 

21. j- det{ A X B] » (detf A X B]) <X* V 

22. jz dei[X'} «idet[X ) • (det[X )) (X‘S # 

23. ^ det f X n ] « n(det(X]) n (X'V 
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